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INTRODUCTION 


The  current  design  trend  toward  higher  tip  speed, 
reduced  engine  weight,  and  higher  pressure  ratios  has  made 
the  modern  high  speed  fan  and  compressor  more  susceptible  to 
supersonic  unstalled  flutter.  Unstalled  flutter,  as 
distinguished  from  stall  flutter,  occurs  at  small  angle  of 
attack  and  with  attached  flows.  In  addition,  it  occurs  at  or 
near  the  design  condition.  Flutter  itself  is  a  self-excited, 
unsteady  phenomenon  which  can  occur  when  the  exciting  forces 
just  equal  the  elastic  forces  of  the  turbomachinery  blades. 
For  a  review  of  the  problem  of  fan  and  compressor  flutter, 
refer  to  Ref.  6. 

The  original  work  in  flutter  analysis  evolved  from  the 
flutter  of  aircraft  wings,  which  was  observed  as  early  as 
1914.  Theodorsen  [Ref.  22]  developed  the  most  basic  analytic 
treatment  of  unsteady,  subsonic  aerodynamics  for  an  airfoil 
oscillating  in  pitch  and  plunge.  This  analysis  dealt  with 
what  has  become  known  as  classical  flutter,  and  showed  that 
flutter  was  only  possible  when  the  natural  frequencies  of 
flexure  and  torsion  coalesced  as  the  relative  velocity  over 
the  wing  was  increased.  Each  mode  of  oscillation  by  itself 
(pitch  or  plunge)  was  stable. 

Cascade  flutter  differs  from  classical  flutter  in  that 
the  inertia  and  stiffness  properties  of  a  turbomachine  blade 
differ  greatly  from  those  of  a  wing.  As  a  result,  the 
critical  flutter  speed,  the  speed  at  which  the  onset  of 
flutter  occurs,  is  very  much  greater.  Garrick  and  Rubinow 
[Ref.  9]  extended  the  work  of  Theodorsen  to  include 
supersonic  flow,  and  their  analysis  predicted  that  single 
degree  of  freedom  torsional  flutter  was  possible,  although 
the  bending  mode  was  still   stable.   Lane   [Ref.  14]   showed 
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that  cascade  flutter  could  be  analyzed  by  considering  only  a 
single  oscillating  blade. 

Work  on  cascades  with  subsonic  leading  edge  locus  is 
fairly  recent.  The  first  attempt  was  made  by  Gorelov 
[Ref.  10].  Verdon  [Ref.  23]  obtained  some  of  the  first 
numerical  results  using  a  finite  difference  procedure.  Bri'x 
and  Platzer  [Refs.  3,4]  were  also  able  to  obtain  results 
using  a  method  of  characteristics  approach.  Another 
solution  was  presented  by  Nagashima  and  Whitehead  [Ref.  15] 
using  dipole  distributions.  All  three  solutions  were  in 
close  agreement.  Sisto  and  Ni  [Ref.  20]  developed  a  'time 
marching"  technique  applicable  to  an  infinite  cascade.  Their 
results  appear  to  be  in  good  agreement  with  the  finite 
cascade  results. 

An  analytic  development,  valid  for  low  frequency  blade 
motion,  has  been  developed  by  Kurosaka  [Ref.  11]  and  is 
currently  being  extended  to  the  higher  frequency  range 
[Ref.  12].  These  results  are  applied  to  an  infinite  cascade. 
Very  recently,  another  analytic  treatment  has  been  presented 
by  Verdon  [Ref.  25]  (see  also  Verdon  and  McCune  [Ref.  24]) 
in  which  an  infinite  series  representation  is  used  to 
express  the  influence  of  neighboring  blades  and  wakes  on  the 
reference  passage.  However,  series  convergence  was  not 
always  obtained,  and  some  doubt  exists  as  to  its  uniform 
validity  [Ref.  13].  A  new  analysis  currently  being 
developed  by  Verdon,  applicable  to  the  infinite  cascade, 
does  appear  to  yield  uniformly  good  results  [Ref.  26], 
Actual  wind  tunnel  tests  of  cascade  flutter  have  been  made 
by   Snyder  and  Commerford  [Ref.  19]  and  by  Fleeter  [Ref.  7]. 

In  the  current  study,  a  linearized  method  of 
characteristics  procedure,  based  on  a  developeraent  by  Teipel 
[Ref.  21]  for  unsteady  aerodynamics  of  a  flat  plat< 
airfoil,   is   used  to  examine  the  pressure  distributions  and 
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flutter  boundaries  for  cascade  blades  having  subsonic 
leading  edge  locus.  It  is  an  extension  of  previous  work  done 
by  Chalkley  [Eef.  5]  and  Brix  [Ref.  2]. 

In  addition,  the  method  of  characteristics  was  applied 
to  the  problem  of  supersonic  flow  through  oscillating 
cylindrical  shells  to  determine  the  pressure  distribution 
and  associated  generalized  forces  acting  on  the  cylinder 
walls. 

Although  it  is  unreasonable  to  suppose  that  the 
two-dimensional  cascade  investigation  can  accurately  model 
the  complex  mixed  transonic-supersonic  flow  of  an  actual 
present  day  turbomachine,  it  is  hoped  that  this  analysis 
will  provide  a  better  understanding  of  the  problems 
involved,  and  serve  as  a  starting  point  for  the  flutter 
analysis  during  the  design  stage  of  modern  high  speed  fans 
and  compressors. 
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A.   PROBLEM  FORMULATION 

Consider  a  cascade  composed  of  a  semi-infinite  array  of 
thin,  two-dimensional  flat  plate  airfoils  immersed  in  the 
supersonic  flow  of  an  inviscid,  non-conducting,  ideal  gas 
having  constant  specific  heats.  Additionally,  the  flow 
field  is  assumed  to  be  irrotational  and  isentropic. 


1  / 

/ 

/  / 

/ 
/ 

>'' 

/  / 

K 

? 

d 

Figure  1.  Cascade  geometry 

The  only  restriction  on  the  cascade  is  that  the  locus  of 
blade  leading  edges  must  be  subsonic;  i.e., 


> a 

2 


(II-  1) 
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where   /?  is  the  stagger  angle  and   a   is   the   Mach   angle 

1 

(  sin"*1—  ).   As  shown  in  figure  1,  the  cascade  is  aligned  in 
H 

the  X-Y  coordinate  system  so  that  the   leading-edge   of   the 

first  blade  is  at  the  origin  and  the  blade  itself  lies  along 

the  X-axis. 


Perturbations  in  the  flow  field  are  caused  by  blade 
motion  in  pitch  and  plunge,  assumed  to  be  of  small  amplitude 
with  simple  harmonic  motion.  For  a  cascade  with  subsonic 
leading-edge  locus,  the  blades  are  swept  aft  of  the  initial 
Mach  wave.  Thus,  disturbances  are  allowed  to  propagate 
throughout  the  domain  of  influence  (the  region  downstream  of 
the  Mach  cone  emanating  from  the  leading  edge  of  the  first 
blade) ,  and  are  not  restricted  to  the  region  between  blade 
passages. 


B.   GOVERNING  EQUATIONS 

Subject  to  the  initial  assumptions  stated  above,  the 
governing  eguations  for  this  boundary  value  problem  are  the 
continuity  eguation 

Dp  <?u    dv 

—   +  Pi—  *   T-)  s  0  (II-  2) 

DT      dX    dY 

the   eguation   of   motion    (   Euler  equations   ) 

Du         1<3p 

(II-3a) 


DT 

+ 

pdx 

" 

0 

Dv 
DT 

+ 

1*2. 

pdY 

= 

0 

(II-3b) 
the  condition  of  irrotationality 
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d\i        dv 

=0  (II-    4) 

ay      ax 

and  the    energy   eguation 

Ds 

™  =  °  (II~  5) 

DT 

D 

—    denotes  the  material,  or  substantive,  derivative   of   a 

DT 

fluid   particle  which  is  a  function  of  position  and  time.  It 

is  composed  of  the  local  derivative  (rate  of  change  with 
respect  to  time)  and  the  convective  derivative  (rate  of 
change  with  respect  to  position)  and  is  written  as 

D(  J    d(   ) 

-zr-  =  "^r~  +  <v°     >  <  >  (n-5a) 

DT    aT 


where  V  is  the  velocity  vector.   The  local  sonic  velocity  is 
given  by, 

aP     dp 
op  s        ap 

and 

P 

— y  =  constant  (II-  7) 

pY 

Taking  the  total  differential  of  Eg.  (7) 

(j)    dp  -  yp(-)      d/>  =  0  (II-  8) 


Rearranging 


dp  p 

—  =y-  (ii-  9) 

dyo  p 
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c2  =r-  (11-10) 

p 


BOUNDARY  CONDITIONS 


The  flow  boundary  condition  (flow  tangency)  for  a   blade 
having  an  equation  of  the  form 

F(X,Y,T)  =  0  (11-11) 

is  developed  as  follows.  At  each  point  of  a  solid-fluid 
surface,  at  every  instant,  the  normal  component  of  the 
relative  velocity  between  the  solid  and  the  fluid  must 
vanish.  In  ether  words,  the  total  time  rate  of  change  in 
F(X,Y,T)  following  a  surface  particle  of  the  solid  is  zero. 
Mathematically,  this  statement  is  expressed  as 

DF 

—  =0  (II-12) 

DT 

For  the  cascade,  the  equations  for  the  upper  and  lower 
surfaces  of  a  blade  are  written 

F  (X,Y,T)  =  Y  -  Y  (X,T)  =  0  (H-13) 


and 

F  (X,Y,T)  =  Y  -  Y  (X,T)  =  0  (11-14) 

L  L 

Applying    Eg.     (12)     to    Egs.     (13)    and    (14) 

DFU       dl0       _dFu         dFy 

u  =   — °+    u— u+    v— u=    0  (11-15) 

DT         dT  <3X  dY 


dT  dX 


(11-16) 


22 


since   —  =  1.   Likewise, 
dl 

df.      ay,    dY{ 

— L  =  v •  -  u— L=  0  (11-17) 

DT        dT  ax 

The  other  boundary  condition  that  must  be  enforced  is  the 
Sommerfeld  radiation  condition,  i.e.,  disturbances  must 
propagate  away  from  their  sources. 


LINEABIZATION 


The  normal  velocity  perturbation  is  written 

v(X,Y  (X,T))  =  — u  +  u— u  (11-18) 

u  dT       ax 

ay,  _ay 

v(X,Y  (X,T))  =  r— U+  u—L  (11-19) 

l  aT       ax 

From  the  assumption  of  thin  blades  and  small  amplitude 
oscillations,  all  flow  quantities  may  be  considered  to  be 
small  quantities  which  are  linearly  superimposed  onto  the 
freestream  quantities.   Thus 

u  =  u  +  u  (11-20) 


c  =  c  +  c  (11-21) 

o 

v  =  v  (11-22) 

while  the  density  and  pressure  perturbations  are 

Ap  =  p  -  p  (11-23) 
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Ap    =    p   -    p  (11-24) 

Also,  because  of  the  assumption  of  thin  blades,   Y    and   Y 

u        L 

are  small   quantities  with  respect  to  axial  distances  along 

the  blade.  Thus,  substituting  Eqs.  (20)  and  (22)  into  Eqs. 
(18)  and  (19)  and  neglecting  the  effect  of  higher  order 
terms,  the  normal  velocity  perturbation  equation  becomes 

v(X,Y  (X,T))  =  —  +  u  — u  (11-25) 

u  dT  odx 


a  y,        a  y, 

V(X#YL(X,I))  -  -♦  u— L  ,11-26, 


To   transfer   the   boundary   condition  to  the  axis,  Egs. 
(25)  and  (26)  are  expanded  in  a  Taylor  series  about   Y  =   0. 

<3v  (X,Y=0+,T) 
V(X,Y  (X,T))  -  v(X,Y  =  0  +  ,T)  +  Y 


dY 

1  a2v(X,Y=0+,T) 

_Y2 .  +  •  •  •        (11-27) 

2  u      (dY)  s 


dv  (X,Y  =  0~,T) 
v(X,Y  (X,T))  =  v(X,Y=0-,I)  +  Y 


L  L       dY 

1  a2v(X,Y=0~,T) 

-Y2  +  •  •  •        (11-28) 

2  L      (dY)2 

Neglecting  higher  order  terms, 

v(X,Y  (X,T))  =  v(X,Y  =  0  +  ,T)  (11-29) 


V(X,Y  (X,T))  =  v(X,Y=0",T)  (11-30) 


and  thus 


dYu    dYu 
v(X,0+,T)  =  — ■♦  u  ■  (H-31) 

dT        oax 


v(X,0~,T)  =  — -L+  u  -— L  (11-32) 


th 
For  cascaded  blades,  it  is  assumed  that   the   n     blade 

st 
leads   the   motion   of  the  (n-1)    blade  by  an  amount  8,  the 

interblade  phase  angle.  Expressing  this  and  the  assumption 
of  sinusoidal  motion,  the  equation  for  the  upper  and  lower 
blade  surfaces  becomes 

i[wT  +  (n-1)Sl 
Y  (X,T)  =  Y  (X) e  (11-33) 


i[wT  +  (n-1)8  J 
Y  (X#T)  =  Y  (X) e  (11-34) 

which   reduces   the   flow   tangency  condition,  Egs.  (31)  and 
(32),  to 

dY  i{(n-1)8}  iwT 

v(X,0+,T)  =  [icuY  (X)  +  u Je         e        (11-35) 

u       °dX 


dY      i{(n-1)8}  iwT 

v(X,0-,T)  =  [icuY  (X)  +  u  ]e         e        (11-36) 

L        °d  X 

iwT 
Now   let    Y  (X,T)   =   Y  (X,T)   =   Y  (X) e    .    Then  the  flow 
u  L 


tangency  condition  for  pitch  is, 


*(X)       =  -  a  (X  -  x  c)  (11-37) 

pitch       o       o 


where, 
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iwT 


(11-38) 


while  for  the  plunge, 

Y(X) 
where, 


plunge 


h  =  h  e 


iwT 


(11-39) 


(11-40) 


Deflections   are   defined  as  positive  downward,  and  angle  of 

twist  positive  leading  edge  up.   h   and  a   are  dimensionless 

o       o 

complex  amplitudes,  with  x   the  twist  axis  location  (elastic 

axis)  in  percent  chord  S.  Figure  2  shows  the  positive 
direction  of  all  quantities  used  in  defining  the  flow 
tangency  condition. 

y 
i 


Figure  2.  Positive  directions  of  Lift 
and  Moment  resulting  from  deflections 
in  pitch  (a)  and/or  plunge  (h)  . 


To   non-dimensionalize   the   flow   tangency    equations, 
distances   are  scaled  to  the  blade  chord  and  time  to  blade 
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chord  divided  by  freestream  velocity  u  .    This  leads  to  the 

o 

introduction  of  the  reduced  frequency  parameter 

CJC 

k  =  (11-41) 

Uo 

which  is  twice  the  reduced  frequency  as  defined  by  Garrick 
and  Rubinow  [Kef.  9]  who  scaled  all  quantities  to  blade 
semi-chord. 

Non-dimensionalizing  Eqs.  (35)  and  (36)   and   using   the 
definition   of   reduced   frequency   above,  the  flow  tangency 
conditions  become 
pitch: 

v(x,0±#t)  ikt 

=   -   a    {cos[  (n-1)  SJ  -    k  (x-x    )  sin[  (n-1)  S  ]}  e 

u  0  0 

0 

ikt 
-   ia    {ksin[  (n-1)  S  J    +    k(x-x  )  cos[  (n-1)  8  ]}  e 

(11-42) 

plunge: 

v(x,0±,t)  ikt 

=    h    {ksin[  (n-1)  S]   -    ikcos[  (n-1)  S  ]}  e 

Uo         ° 

(11-43) 

To  linearize  the  expression  for  the  pressure,  Eqs.  (10)  , 
(20)  and  (24)  are  used.  Expand  Eg.  (10)  in  terms  of 
perturbation  quantities.  Through  Eg.  (48) ,  second  order 
terms  are  neglected  whenever  they  occur. 

(p  +  Ap) 
cz  +  2cc   =  y— s K —  (11-44) 

P 

o 

Ap 
As  — —  <  1  ,  the  denominator  can  be  expanded  in  a  geometric 
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series    of   the   form 


1 

=1-      c+€2-€3+#e«  (11-45) 

1      +     € 


to    yield, 


c2    +    2cc      =    y(p      +     Ap)  (1 £)  (11-46) 

o  o         '    ro  r  p 

Expanding, 

c2  Ap 
c2    +    2cc      =       (1 e-— -)Ap  (11-47) 

0  0  y  Ap 

However,    as     Ap   and    Ap    are    both   small   quantities 

Ap  dp  -1  1 

-£  =    (-)         =    —  .  (11-48) 

Ap  dp  c2 

Expanding   Eg.   (47)   ,   substituting   into   Eq.   (46)    and 
rearranging  terms, 

Ap     2  2yc        -» 

=  cc  [  1  +  l  (11-49) 

Po    r-l   oL     (X-1)c0J 

2yC 

But  <   1.    Using   Eg.   (43)  to  expand  the  term  in 

(r-Dco 

brackets   yields   the   linearized   form   of    the    pressure 
difference 

2 

p  -  p   =  p  c  (c  -  c  )  (11-50) 

r  0     y-  1  r0  0         0 
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III.   METHOD  OF  CHARACTERISTICS 

A.   CHARACTERISTIC  DIRECTIONS 

To  introduce  the  method  of  characteristics,  it  is 
important  to  recall  that  for  supersonic  flow  there  are 
distinct  regions  of  influence  and  dependence.  Hence,  curves 
must  exist  along  which  the  values  of  the  flow  variables  and 
their  derivatives  are  known,  yet  the  flow  field  in  the 
immediate  neighborhood  of  these  curves  is  indeterminate. 
Therefore,  to  obtain  the  directions  in  the  x-y  plane  of 
these  curves,  which  are  called  characteristic  directions,  or 
more  simply,  characteristics,  one  must  prescribe  that  the 
derivatives  across  these  curves  are  indeterminate. 

Mathematically,  this  means  that  the  dependent  variables 
must  be  continuous,  but  their  first,  or  higher  order, 
derivatives  may  be  discontinuous.  Physically,  the 
characteristics  act  as  information  carriers,  in  that  ail 
disturbances  propagate  along  them.  In  particular,  they 
represent  the  paths  of  Mach  waves,  and  are  therefore 
referred  to  as  Mach  lines. 

To  facilitate  the  analysis,  it  is  desirable  to  introduce 
new  coordinates  £  and  7)  .  The  governing  equations  are 
written  in  terms  of  the  two  arbitrary,  but  intersecting 
curves 

£  =  £  (X/Y)  =  constant  (III-1a) 

and 

V   =  ?!  (x,y)  =  constant  (III-1b) 

However,  the  equations  are  more  suitably  transformed  to  the 
arbitrary  coordinates  if  the  dependent  variables   are   u,   v 


29 


and      c".         Here,      air      is     assumed      to      be   a    perfect   gas   with 

constant      specific      heats      c         (constant      pressure)       and      c 

p  v 

(constant  volume).  Rewriting  Eg.  (11-10)  as 


yRQ  (III-  2) 


R  =  c   -  c  (III-  3) 

p     v 

where    R    is  the   gas   constant   for   air,   and    0    is 
temperature   in    degrees    Rankine.     Taking    the   total 

differential  of  Eg.  (2)  and  dividing  by  c2  , 

dc    de 
2~r  =  —  (III-  4) 

c    e 

Taking  the  total  differential  of  Eg.  (11-10), 

1       do 
2cdc  =  v-dp  -   p—  (III-  5) 

or,  dividing  by  c2, 

dc    dp   dp 
2—  =  —  ~  —  (III-  6) 

c    P     p 

Now  from  the  First  Law  of  Thermodynamics  and  the   fact   that 
the  flow  is  homentropic, 

p 

Gds  =  du dp  (III-  7) 

int    pz   r 

For   an   ideal  gas,  du    =  c  dG   and  Eg.  (7)  may  be  written 
int     v 


dp 
Gds  =  c  de  -  R9—  (III-  8) 

v        p 
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Dividing  by  c  9 
v 

1       de        dp 

—  ds  = (y-1)—  (III-  9) 

cv     e        p 
and  using  Egs.  (4)  and  (6)  yields 

1        dc        1 

—  ds  =  2 (/-1)-dp  (111-10) 

cv       c         p 

But  because  the  thermodynamic  properties  of  a  continuous 
flow  field  are  to  be  considered  as  continuous,  single- valued 
scalar  point  functions,  Eg.  (10)  must  hold  for  any  change  of 

the  fluid  properties  s,  c  and  p.  Hence,  it  may  be  written  in 
terms  of  material  derivatives.  Then,  subject  to  Eg.  (II-5)  , 
Eg.  (10)  becomes 

1  Ds      1  Dc         1Dp 

=  2  — (r-1) =  0  (111-11) 

C  DT      C  DT         pDT 

v  c 

Then,  using  Eg.  (11-10) ,  Eg.  (11)  reduces  to 

Dp     2   Dc 

—  =  c—  (111-12) 

DT    y- 1   DT 

Now  the  perturbation  eguations  are  introduced.  Writing  Eg. 
(II-9)  as 

Dp    _  Dp 

=  cs—  (111-13) 

DT      DT 

and  substituting  into  Eg.  (12)  allows  Eg.  (II-2)  to  be 
expressed  as 

2  dc  2        dc  du  dv 

+ ,u  +    c +    c  =    0  (111-14) 

y-1dT       7-1    °dx  odx  <>dY 

neglecting  the  effect  of  higher  order  terms.  This  is  the 
form  of  the  continuity  eguation  which  is  to  be  used  in  the 
transformation. 
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To  render  the  first  Euler  equation,  Eq.  (II-3) ,  in  a 
form  suitable  for  transformation,  use  Eg.  (11)  in 
differential  form 

2-^-dc  =  (X-1)-dp  (I  II- 15) 

c  p 

With  Eg.  (11-9)  this  yields 

1       2 

-dp  =  cdc  (111-16) 

p         y-i 

Expanding  the  total  derivatives  and  equating  the  partial 
derivatives,  one  obtains 

1aP    2  dc 

=  c" —  (111-17) 

Pdx     x-i  ax 

Substituting  into  the  first  Euler  eguation 

Du"    du   _du    du  2   _dc 

—  =  —  +  u  —  *   v—-  = c  —  (111-18) 

dt  3t        ax    ay     r-1  ax 

Introducing  the  small  perturbation  guantities  and  neglecting 
second  order  guantities  yields  the  desired  form  of  the  first 
Euler  eguation, 

an    an   2   ac 

__  +  u  __  +  c  __  =  0  (111-19) 

aT       odx      y--[  oax 


The  condition  of  irrotationality 

av      an 

-  -  -  =  0  dii-20) 

is  used  in  place  of  the  second  Euler  equation  as  one  of  the 
governing  equations  to  be  transformed.  Then  the  continuity 
equation,  the  first  Euler  equation  and  the  condition  of 
irrotationality  form  a  system   of   three   equations   in   the 

three   unknowns   u,   v  and   c.  For  the  assumption  of  simple 


3  2 


harmonic  motion,  the  Teipei  amplitude   functions   [Ref.   21] 
are  introduced, 


U(X,Y)e 


rwT    u  -  u. 


(III-21a) 


V(X,Y)  e 


(III-21b) 


C(X,Y)e 


icuT 


2   1c- 

H  ?   c~ 


(III-21c) 


where  U,  V  and  C  are  complex  non-dimensional  amplitudes.  U 
and  V  are  the  perturbation  velocities  and  C  is  the 
perturbation  in  sonic  velocity  which  will  be  related  to  the 
pressure  coefficient  through  Eg.  (11-49) . 


Substituting   Egs.   (21)  into  Egs.   (14),  (19)  and  (20), 
and  then  non-dimensionalizing  as  before  yields. 


au         t dV  dc 

—  +    /M2-1  —  +   »2 —  +    ikM2c 
dx  dy  dx 


(III-22a) 


au      dC 

—  +   —  +   ikU   =   0 
dx        5x 


(III-22b) 


au       , av 

/M2--|    =o 

ay  dx 


(III-22c) 


For     £    =   constant   or    77     =   constant,    both     £  (x,y(x))      and 
i?(x#y(x))    are   implicit   functions   of      x.      Hence, 


d£    =  £    dx   +  £    dy   =   0 
x  y 


(III-23a) 


and 
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drj   =   7) 
X 

dx 

+ 

y 

=   0 

Therefore, 

X__          X 

ey      ^y 

= 

- 

_dy 
dx 

(III-23b) 


(III-  24) 


and  the  labeling  of  the  characteristic  directions  is 
arbitrary.  By  convention  the  left-running  Mach  line  is 
normally  denoted  £   and  the  right-running  Mach  line,  77  • 


To  solve  for  the  characteristic  directions  in  the  x-y 
plane,  Egs.  (22)  are  rewritten  along  lines  of  constant  £ 
and  77. 


£   U      +    JK^A   £    V      +    M2£   C 


ikM2C 


y  i 


(III-25a) 


C    u      +    C  c 


u      - 


ikU 


x   75 


(III-25b) 


(     U         -      /M2-1      {     V         =     -     77     U         +      /M2-1      77     V 

y  x  y  77  x  77 


(III-25c) 


Now   Cramer's    rule    is   applied   to   solve   for   the 

derivatives  of  the  dependent  variables   with   respect   to   £ 

along   77  =   constant.   To   satisfy  the  condition  that  these 

0 
derivatives  be  finite  requires  that  each  be  of  the   form  — , 

0 

i.e.,         indeterminate.         Equating      the      determinant      of      the 

t 
coefficient    matrix    of    {   U    ,    C    ,    V      }       to   zero   yields, 

£      E   £2(K2-1)    "  £2    1    =    0  (III-    26) 

xx  y 
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to  which  there  are  three  real  solutions,  each  of  which 
defines  a  characteristic  direction  in  the  x-^y  plane.  These 
are, 


dy 
dx 


e   =  o 


(III-27a) 


dy 

dx 


1 


(III-27b) 


The   same  results  are  obtained  when  one  solves  for  the  first 
derivatives  with  respect  to  77  along  lines  of  £  =  constant. 


Eg.   (27a)  gives  the  streamline  slope  in  the  flow  field, 
while  Eg.  (27b)  gives  the  slope   of   the   two  curves   £ 
constant   and   17   =  constant.  It  should  be  noted  that  these 
are  the   same   characteristic   directions   as   obtained   for 
steady,  two-dimensional  flow.   By  convention, 


dy  1 

dx  £  =const      /M2-1 


(III-28a) 


dy 

dx  i7=const 


1 


(III-28b) 


dy 

dx  stream 


(III-28c) 


For    an   arbitrary   function   f   =   f  (x,y)  ,   the   total 
derivative  is 


df     df 

df  =  — dx  * dy 

dx     dy 


(III-  29) 


For  f(x,y)  taken  along  a  curve  as  defined  by  Egs.  (1), 
however,  y  is  an  implicit  function  of  x  and  Eg.  (29)  is 
rewritten  as 


df  _  df        df dy 
dx   dx   dydx 


(III-  30) 


Therefore,  along  a  streamline, 


df 


df 


dx  stream    dx 


(III-  31) 


while  along  the  Mach  lines. 


df 


d£      1  df 


dx  £  =const  dx      /H2-1    dy 


(III-32a) 


df 
(— ) 


df 


1      df 


dx77=const  dx      /M2-1    dy 


(III-32b) 


Then,   solving   for   the  partial  derivatives  in  terms  of  the 
ordinary  derivatives, 


df    df 

dx    dx  str 


(III-33a) 


df    1   df  df 

dx    2   dx  £  =const    dx  rj  =const 


(III-33b) 


df         1   , df  df 

— ■    =   — /H2-1     [  (— )  -      ( )  ] 

dy        2  dx  £  =const  dx  7j=const 


(III-33c) 


The  compatibility  relations  are  obtained  by  writing  Egs. 
(22)  using  Egs.  (33) .  Thus, 


1  dU      dU       1         dV 

2  dx  £     dx  77     2         dx  £ 

H2   dC       dC 


dV 
ax  77 


( — )  ]  +  ikM2C  =  0 
dx  77 


(III-34a) 
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1  dO      dU       1   dC      dC 

of  <>Th    +  <7">  3  +  ~C  (— ) .  +  (t- ■)  ]  +  ikp  =  o 

2  dx(     dx  7)  2   dx  4     dx  rj 

(III-34b) 

dU      dU        dV       dV 
dx  £     dx  -17       dx  £     dx  7; 

Eg.  (34b)  is  multiplied  by  M2  and  then  subtracted  from  Eg. 
(34a) .  The  resulting  eguation  is  first  added  to,  and  then 
subtracted  from  Eg.  (34c)  to  give  two  of  the  compatibility 
equations 

dU       dV         M2 

<— )   -  (— )   ♦  ik— —  (U  -  C)  =  0  <III-35a) 

dx  £  dx  £      M2-1 

dU       dV         M2 

( )   +  ( — )  *    ik (0  -  C)  =  0  (III-35b) 

dx  17    dx  7)  M2-1 

The  third  compatibility  relation  is  obtained  from  writing 
Eg.  (34b)  along  a  streamline. 

dU  dC 

( — )       +  ( — )        +  ikU  =  0  (III-35c) 

dx  stream     dx  stream 

The  governing  eguations  have  now  been  reduced  to  a  system  of 
three  ordinary  differential  equations  which  are  implicitly  a 
function  of  x,  with  U  ,  V  and  C  complex  amplitudes.  Thus, 
Eqs.  (35)  are  really  six  eguations,  with  three  real  and 
three  imaginary  components. 

To   write  the  boundary  conditions  in  terms  of  the  Teipel 
amplitude  functions,  note  that 

. ikt 

v(x,0,t)  =  u  yP^T  V(x,y)e  (III-  36) 

from  the  definition  of  the  amplitude  function,  and  therefore 
the  flow  tangency  condition  is, 
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V  =  V  (x,y)  = 


1   v(x,0,t) 
'M2-1     U„ 


(III-  37) 


v  (x,0,t) 
where  is  given  by  (11-40)  for  pitch  and  (11-4  1)  for 

uo 

plunge. 


B.   FINITE  DIFFERENCE  EQUATIONS 

To  write  Egs.  (111-34)  in  finite   difference   form,   the 
computational  molecule  shown  in  figure  3  is  used. 


current     Mach 
line 


previous    Mach 
F       \  line 

12 


Figure  3.   Computational  molecule  for  a 
general  flow  field  point. 


If  F  is  an  arbitrary  complex  flow  guantity,  derivatives 
with  respect  to  x  along  the  characteristics,  in  finite 
difference  fcrin,  are  written, 


dF 

dx  str 


'22 


JJ_ 


Ax 


(III-38a) 
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dF       F    -  F 

tTK  =  I —  (III-38b) 

dF       F99  -  ?0. 

( — )  =  -&£- ±±-  (III-38C) 


The  values  of  the  perturbation  quantities  in   Eqs.   (38) 
are  averages,  so  that 

1 
(F)     =  -(  F    ♦  F    )  (III-39a) 

str    2    11     22 


1 
(F)   =  ■  — (  F    *  F    )  (III~39b) 

£    2    12     22 


(F)   =  o(  ?oo    +  Foi  }  (III-39C) 

rj  2         22  21 


1 .   Gen eral  Flow  Field  Equations 

With  these  relations  substituted  into  Eqs.  (35) ,  the   system 
of  equations  becomes, 

U     +  C    -AU     =K  (III-40a) 

22R     22R     I  221     12R 


0     +  C     +AU     =K  (III-U0b) 

221     221     I  22R     121 


U     -V     -  B  (  U     -C     )=K  (III-40C) 

22R     22R     I    221     221       34R 


U     -V     +  B  (  U     -  C     )=K  (III-aOd) 

221     221     I    22R     22R       341 


U     +  V     -  B  (  0     -C     )=K  (III-40e) 

22R     22R     I    221     221       56R 


39 


U     ♦  V     +B(U     -  C     )=K 
221     221     I    22R     22R       561 


<III-40f) 


where 


K     =U    +  C    ♦  A  U 
12B     11R     11R     I  111 


(III-UOg) 


K     =U     +  C     -  A  U 
121     111     111     I  11R 


(III-4Qh) 


K     =U     -V     +B(U     -  C     ) 
34R     12R     12R     I    121     121 


(III-40i) 


K     =U     -V     -B(U     -  C     ) 
341     121     121     I    12R     12R 


(III-40J) 


K     =U     +V     +B(U     -C     ) 
56R     21R     21R     I    211     211 


(III-40k) 


K     =U     +  V     -B(U     -C     ) 
561     211     211     I    21R     21R 


(III-401) 


and 


A   =  -kAx 
I    2 


(III-41a) 


1    M2 

B   =  -k( )AX  (III-41b) 

I    4   M2-1 

As  shown  by  Teipel  [Ref.  21],  solving  these  equations  gives 


(  1  -ArB.)|A(  K-AD+Kt,ca)-  BTK,OT    I  +2Bf|4(K-.t+K_et)*BTK. 


U 


22R 


221 


ItV[2^34R**56R)~BIK12I  j  ^[gn^  *561  j  +  »ik12rJ 


(Iir-42a) 


(1-W3(K34I  +  K56T)-BIK12rJ 


2B 


J2( K34R*  K5 


(1-A.BT)2  ♦    (2BT)2 


1 

V  =  -(    K  -    K  ) 

22R         2         56R  34R 


(III-42b) 
(III-42C) 


V     =  -(  K     -  K     ) 
221    2    561     341 


(III-42d) 


C     =K     -U     +  A  U 
22R     12R     22R     I  221 


(III-42e) 


C     =K     -  U     -  A  U 
221     121     221     I  22R 


(III-42f) 


These   are   the  finite  difference  equations  for  a  grid  point 
located  in  the  general  flow  field. 

2.   U££<2r  Surface  of  a  Blade 


For  a  grid  point  on  top  of  a  blade,  the 
computational  molecule  of  figure  3  is  modified  as  shown  in 
figure  4. 

y 


Figure  4.   Computational  molecule  for  a  grid 
point  on  the  upper  surface  of  a  blade. 

Perturbation   quantities   at   grid   point  F    are  considered 

zero  and  the  normal  velocity  V    (pitch  or   plunge   at   grid 

point   F   )   is  prescribed  by  the  blade  movement  as  given  by 
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Eg.   (11-37) .   Here  the  applicable  equations  are, 


U     +C    -AU     =K 
22R     22R     I  221     12R 


U     ♦  C     +AU     =K 
221     221     I  22R     121 


U     -BU     +BC     =  K     -  K 
22R     I  221     I  221     56R     34R 


U     +BU     -  B  C     =K     -K 
221     I  22R     I  22R     561     341 


(III-43a) 
(III-43b) 
(III-43c) 
(III-43d) 


where   K    and   K    are  as  before,  but 
12        56 


K     =  V 
34R     22R 


(III-43e) 


K     =  V 
341     221 


(III-43f) 


as  given  by  the  flow  tangency  condition.  Blade  position  is 
always  given  relative  to  the  leading  edge  of  each  blade. 
Solving, 


u  _  "-AIBI)1K56R-K34R-B,K1?1J-2B,|K56,  -K,4,'B,K1?rJ 

22R  (1- ATBT)2  *   (2BT)2 


(II  I -44a) 


(1-AjBj)2  +   (ZBj)2 


C     =K     -U     +AU 
22R     12R     22R     I  221 


(III-44b) 
(III-44C) 


C     =  K     -  0     -AU 
221     121     221     I  22R 


(III-44d) 


and  V    are  known. 
22 
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3.   Lower  Surface  of  a  Blade 


At  a  grid  point  on  the  lower  surface  of  a  blade,  the 

normal  perturbation  velocity   at   grid   point   F    is   also 

22 

known.    The   computational  molecule  used  is  shown  in  figure 
5. 


Figure  5.   Computational  molecule  used  for  a 
grid  point  on  the  lower  surface  of  a  blade. 


For  this  grid  point,  the  finite  difference  equations 


become, 


0     +C     -  A  U     =  K 
22R  ,   22R     I  221     12E 


(III-45a) 


U     +C     ♦  A  U     =K 
221     221     I  22R     121 


(III-45b) 


U     -  B  (  U     -  C     )=K     +  K 
22B     I    221     221       56R     34R 


(III-45c) 


U     +  B  (  U     -C     )=K     +K 
221     I    22R     22R       561     341 


(III-45d) 


where   K     is  the  same  as  Egs.  (39a)  and  (39b)  ,  but 
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K     ?  0     -  V  *   +  B  (U     -C    ) 
56R     12R     12R     I   121     121 


(III-46a) 


K  =U  -V  -    B     (U  -    C         ) 

561  121  121  I       12R  12R 


(III-46b) 


K  =    V 

34R  22R 


(III-U6c) 


K  =    V 

341  221 


(III-46d) 


Solving   for    the   perturbation    quantities      at      grid      point      22 
yields, 


(1-AlBI)LK56R"K34R-BIK12lJ^2Blb5 


llK56I  *  K3 


221 


!JWJ 


(1-A.BJ     KRCT  ♦  K,„T  ♦  BTK™  -2B,  KKCD*  K,-D-  BTK,. 


561        341  T  "P12RJ         IL  56R       34R       I    12lJ 
(1  -AjBj)2  +    {ZBJ* 


C  =    K  -    0  +AU 

22R  12R  22R  I    221 


(III-47a) 

(III-47b) 
(III-U7c) 


C  =    K  -    0  -    A    U 

221  121  221  I    22R 


(III-47d) 


and   V        are    known. 
22 


4 •      Initial   Left   Running    Mach   Line 


To  obtain  the  equations  for  the  initial  left  running 
Mach  line  (£  =  constant),  the  computational  molecule  of 
figure    5    is    used.    However,    it    is    assumed      that      grid      points 

F  and      F  are      just      in    the   freestream   and   the    limiting 

11  21 

condition      x-^0   is    taken.    For      x— 0,    A      and   B   -»-0.    Thus   Egs. 

I  I 


(38)  reduce  to, 


U    ♦  V    =0  (III>48a) 

22     22 


U    +  C    =0  (III-48b) 

22     22 

or, 

U  =  -  v  =  -  C  (III-  49) 

Now  Eg.  (35a)  becomes 

dU  M2 

(  — )   =  -  ik u  (III-  50) 

dx  £  M2-1 

which  can  be  integrated  along  £  =  constant.  Or, 

U  =  Ul    exp[-ik x]  (III'  51) 

lx=0       M2-1 

Ul     is  known  from  Egs.  (49)  and  (37)  Thus, 
|x=0 

M2  M2 

U    =  -  V    (O)cos(k x)  -  V    (O)sin(k x) 

22R       22R         M2-1       221         M*-1 


(III-52a) 


U    =  V    (O)sin(k x)  -  V    (O)cos(k x) 

221     22R         M2-1       221         M*-1 

(III-52b) 

with  the  appropriate  form  of  Eg.  (37)  used  for  V   (0).  Then, 

22 

V     =  -  U  (III-52C) 

22R       22R 


V     =  -  U  (III-52d) 

221       221 


C      =.  -  0  (III-520) 

22R       22R 
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C     =  -  U  (III-52f) 

221       221 

Eqs.  (52)  are  also  used  to  compute  the  perturbation 
quantities  at  the  leading  edge  of  the  upper  surface  of  the 
first  blade. 

5-   Initial  Right  Running  Mach  Line 

The  developement   for   the   perturbation   quantities 

along   the  initial  right  running  Mach  line  (17=  constant)  is 

analagous  to  the  development  for  the  left.  Here,  F    and  F 

12       11 

of   figure   2  are  assumed  to  be  just  in  the  freestream.  Eqs. 

(38)  reduce  to, 

U   -  V   =0  (III-53a) 

22     22 


U    +  C    =0  (III-53b) 

22     22 


U  =  V  =  -  C  (III-  54) 

Then  Eg.  (35b)  becomes, 

dU  M2 

( — )   =  -  ik x  (III-  55) 

dx  77        M2-1 

Integrating  along  77  =  constant  yields, 

M2 

U  =  Ul    exp[-ik x]  (III-  56) 

|x=0    L    M2-1 


M2  M2 

U    =  V    (O)cos(k x)  +  V    (O)sin(k x) 

22R     22R         M2-1       221         «2-1 

(III-57a) 
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U     =  -  V    (O)sin(k x)  +  V    (O)cos(k x) 

221       22R         M2-1       221         M2-1 

(III-57b) 

Again,  the  appropriate  form  of  Eq.  (37)  is   used   to   define 

V   (0)  .  Then, 
22 

V     =  U  (III-57c) 

22R     22R 


V  =0  (III-57d) 

221     221 


C     =  -  U  (III-57e) 

22R       22R 


C     =  -  U  (III-57f) 

221       221 

Egs.  (57)  are  used  to  compute  the  perturbation  quantities  at 
the  first  grid  point  on  the  lower  surface  of  the  first 
blade.  In  the  computer  program,  lower  surface  points  at 
grid  point  F   are  subscripted  33    or  55  for  pitch  and  plunge 

respectively. 

6.   First  Grid  Point  on  a  New  Blade 


For  the  first  grid  point  encountered  on   any   blade 

other  than  the  first,  the  computational  molecule  of  figure  6 

is   used.    The   new   blade   is   assumed   to   protrude   an 

inf initisimal   distance  past  grid  point  F    which  is  labeled 

22 

F     and  F    in  the  figure.  Then,  the   equations   developed 
22U       22L 

for   grid   points   on  the  upper  and  lower  surface  of  a  blade 

are  used  to  compute  the  perturbation  quantities  at  F     and 
.  22U 


F    respectively. 
22L 
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Figure  6.   Computational  molecule  used  for 
the  first  grid  point  on  a  new  blade. 


(7)  Wake  Grid  Points 

y 


F2. 

blade           Fny 

\  F22U 

slip 

F..L 

F,2 

/    F22L 

plane 

Figure  7.   Computational  molecule  used  for 
a  grid  point  in  the  wake. 


It   is   assumed  that  the  first  wake  grid  point  is  an 
inf initismal  distance  aft  of  the   blade.   To   calculate   the 


flow   field   quantities   in   the   wake,   the  computational 

molecule  is  divided  into  two  triangular  molecules  as   shown. 

Then,   the  equations  developed  for  a  grid  point  on  the  upper 

surface  of  a  blade  apply  to  F     while  those  for   the   lower 

22U 

surface   of   a   blade  apply   to  F    .    The  computational 

22L 

molecule  is  shown  in  figure  7. 


U      ♦  C      -  A  U      =  K  (III-58a) 

22RU     22RU     I  22IU     12RU 


U      +C      +  A  U      =K  (III-58b) 

22IU     22IU     I  22RU     12IU 


U      +V      +B(C      -  U      )=K       (III-58C) 
22RU     22RU     I    22IU     22IU       56R 


U  +V      -  B  (C      -  U     )=K         (III-58d) 
22IU     22IU     I   22RU     22RU      561 

and, 

U  +  C      -AU      =K  (III-59a) 

22RL     22RL     I  22IL     12RL 


U      +  C      +  A  U      =  K  (III-59b) 

22IL     22IL     I  22RL     12IL 


U      -V      +B(C      -U      )=K       (III-5SC) 
22RL     22RL     I    22IL     22IL       3UR 


U      -V      -B(C      -  U      )  =  K       (III-59d) 
22IL     22IL     I    22RL     22RL       341 

wher.e, 

K      =U      +C      +AU  (III-60a) 

12RU     11RU     11RU     I  11IU 


K      =U      +C      -AU  (III-60b) 

12IU     11IU     11IU     I  11RU 


K      =U      +  C      +  A  U  (III-60C) 

12RL     11RL     11RL     I  11IL 


K      =  0      +C      -  A  U  (III-60d) 

12IL     11IL     11IL     I  11RL 


K     =U     -V     +B(U     -  C     )        (III-60e) 
34R     12R     12R     I    121     121 


K     =U     -V     -  B  (  U     -  C     )         (III-60f) 
341     121     121     I    12R     12R 


K     =  0     +  V     +B(U     -C     )        (III-60g) 
56R     21R     21R     I    211     211 


K     =  0     +V     -B(U     -C     )         (III^60h) 
561     211     211     I    21R     21R 


Eqs.  (58)  *■  (60)  form  a  system  of  eight  equations 
with  twelve  unknowns.  However,  in  the  linearized  problem, 
the  wake  cannot  support  a  pressure  difference,  and  the 
normal  velocity  across  the  wake  must  be  a  constant.  As  a 
result  of  these  two  conditions,  to  first  order  effect,  there 
will  be  no  discontinuities  reflected  from  the  wake. 
Therefore,  the  relations, 

V      =  V      =  V  (III-61a) 

22RU     22RL     22R 


V      =  V      =  V  (III-61b) 

22IU     22IL     221 


C      =  C      =  C  (III-61C) 

22RU     22RL     22R 


C      =  C      =  C  (III-61d) 

22IU     22IL     221 

are  used  to  reduce  Egs  (58)  and  (59)  to  a  determinant   eight 
by  eight  system  of  equations  which  are: 
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and 


U      +C     -  A  U      =  K  (IIJ>62a) 

22RU     22R     I  22IU     12RU 


U      +C     ♦  A  U      =K  (III-62b) 

22ia     221     I  22RU     12IU 


U      +V     +B(C     -U      )=K         (III-62C) 
22RU     22R     I    221     22IU       56R 


U      +  V     -  B  (  C     -  U      )=K  (III-62d) 

22RU     221     I    22R     22RU       561 


U      +C     -AU      =K  (III-63a) 

22RL     22R     I  22IL     12RL 


U      +C     +  A  U      =  K  (III-63b) 

22IL     221     I  22RL     12IL 


U      +V     +  B  (  C     -U      )=K         (III-63C) 
22RL     22R     I    221     22IL       34R 


U      +V     -B(C     -U      )=K  (III-63d) 

22IL     221     I    22R     22RL       341 

Egs.  (60)  remain  unchanged. 


C.   THE  PRESSURE  COEFFICIENT 

The   non-dimensional   pressure   distribution   along   the 

surface  of   a   blade   is   determined   using  the   linearized 

pressure   difference  of   Eg.   (11-48)   and  Teipel's   sonic 

velocity   perturbation   amplitude   function,  Eg.     (21c) . 
Rewriting  Eg.  (11-48)  as, 

P  -  P   =  -p    c2(- ^)  (III-  64) 

o   y-Y  o  o  c 

and      then    dividing    by    the   freestream    dynamic  pressure    yields 
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p  -  p  ikt 

- °=  2C(x,y)e  (III-  65) 

Hence,  the  complex  pressure  difference   across  a   blade  at 

grid  point  F   resulting  from  oscillation  in  pitch  or  plunge 
22 

is, 

p   -  p  ikt 

P(x,y,t)  =  -^ -u  =  2  (C    -  C    )e  (III-66a) 

|yPM2        22L     22U 


This  expression  for  the  pressure  difference  across  a 
blade  appears  to  be  independent  of  interblade  phase  angle  8. 
However,  the  interblade  phase  angle  was  used  for  the 
computation  of  the  normal  velocity  through  the  flow  tangency 
eguations,  and  they  were  used  in  the   determination   of   C. 

Therefore,   to   compute   the  pressure  distribution  along  the 

th 
n    blade  when  it  is  in  a  specific   position   requires   that 

the   position   be   specified  in  terms  of  the  position  of  the 

first  blade.  In  general,  the  complex  pressure  is 

ikt 
P(x,y,t)  =  (P   +  iP  )  e  (III-  67) 

R      I 

th 
with  P   and  P   defined  by  Eg.  (66)  .  The  position  of  the  n 
R       I 

blade  is 

(kt)   =  (kt)  +  (n-1)8  (III-  68) 

n 

Then  expanding  Eg.  (67)  , 

P(x,y,t)  =  [P  cos(kt)  -  P  sin(kt)  ] 
R  I 

i[P  sin(kt)  ♦  P  cos(kt)  1       (III-  69) 
R  I 
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Figure  8  shows  a  complete  oscillation  cycle  for  a  blade 
undergoing  simple  hamonic  motion  in  pitch.  The  extension  to 
plunge  is  cbvious,  and  the  following  discussion  is 
applicable  there  also. 


kt  = 

0 

kt  = 

90 

kt  = 

180 

kt  =  270 

Figure  8.   Blade  oscillation  cycle. 
The  two  conditions  of  primary  interest  here  are  when  the 


th 
n    blade  is  in  the  initial  position 

(kt)   =  0 


(III-  70) 


or  in  the  mean  position,  pitching  leading  edge  up, 

(kt)   =  270  (III-  71) 

n 

Then  solving  for  the  position  of  the  first   blade   when   the 

th 
n    blade  is  in  the  initial  position  gives, 

(kt)  =  -  (n-1)S  (III-  72) 

and  in  the  mean  position, 

(kt)  =  270  -  (n-1)S  (III-  73) 

Except  for  flutter  analysis,  the  real  component  of  the 
pressure  is  normally  of  primary  interest,  and  can  be  found 
by   solving   Eg.   (69)   with  (kt)  expressed  as  a  function  of 
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(kt)  .   For   ease   of   showing   the   relation   of   the   real 
n 

component   of   pressure   for   the  two  blade  positions  above, 
consider  the  special  case  of  8=0.  Then, 


Re{P(x,y,t)}      =  P 
kt=0     B 


(III-74a) 


Re{P(x,y,t)}       =  P 
kt=270    I 


and  thus, 

Re{P(x,y,t)} 


=  Im{P(x,y,t)} 
kt=0  kt=270 


(III-74b) 


(III-75a) 


Re£P(x,y,t)}        =  Im{P(x,y,t)} 

kt=270  2  kt=0 


(III-75b) 


Having  obtained  the  pressure  distribution  over  the 
blades,  the  dimensionless  lift  per  unit  span  and  the  moment 
per  unit  span  square  are  found  from, 


L(x 


,y,t)  =  =   /P(x,y, 

■Jn 


t)  dx 


and 


M(x,y,t)  =  =   /(x-x  )P(x,y,t)dx 
Jo  ° 


(IIl>76a) 


(III-76a) 


where  L(x,y,t)  and  M(x,y,t)  are  complex  quantities  which 
include  the  effects  of  blade  oscillation  in  pitch  and 
plunge.  Lift  is  positive  upward  and  moments  are  positive 
clockwise  (see  figure  2) . 
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IV.  FLUTTER  ANALYSIS 


A.   EQUATIONS  OF  MOTION  FOR  A  SINGLE  BLADE 

Lane  [Ref.  14]  showed   that   cascade   flutter   could  be 

analyzed   by   considering   only   a  single  oscillating  blade. 

Following  Fung  [Ref.  8],   the  eguations  of  motion  for  a  flat 

plate   airfoil  immersed  in  an  inviscid  flow  are  derived  from 

a  force  balance  of   the   inertia,   elastic   and   aerodynamic 

forces.  For  bending  motion,  h  is  measured  positive  down  from 

the  elastic  axis.  The  torsional  motion  is  positive  clockwise 

about   the   elastic   axis  (refer  to  figure  2) .  It  is  assumed 

that  a  pair  of  springs,  with  spring   constants   K    and  K  , 

a        h 

oppose  the  blade  motion. 

The  total  inertia  force  per  unit  span  is, 

F   =  -  (M  h  +  S  a)  (IV-  1) 

i       b     b 

and  this  force  exerts  a  moment  about  the  elastic  axis  x  of, 

o 

M   =  -  (I  a  +  S  h)  (IV-  2) 

i       a      b 

where 

M   -  total  blade  mass  per  unit  span 
b 

S   -  wing  static  moment  (about  x  ) 
b  o 

I   -  wing  mass  moment  of  inertia  (about  x  ) 
a  o 

The  springs  produce  restoring  forces  and  moments  of 

F   =  -  hK  (IV-3a) 

r       h 
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H   =  -  aK  (IV-3b) 

r       a 

Then  the  equations  of  motion  for  the  blade  are, 

Lift/span  =  M  h  +  S  a  +  hK  (IV-4a) 

b     b      h 


Moment/span  =  S  h  +  I  a  +  aK  (IV-4b) 

b     ft      a 

The  assumption  of  simple  harmonic  motion  is  expressed  as, 

iwT 
h  =  h  e  (IV-5a) 


icuT 
a  =  a  e  (IV-5b) 

o 

where    h    and  a        are   complex   amplitudes.    The   blades 
o         o 

oscillate  with  a  constant  interblade  phase  angle  8 .  A  phase 

difference  exists  between  twist  and  flexure  also,  but  cannot 
be  determined  until  a  flutter  solution  is  obtained.  Then, 
the  ratio  of  blade  deflection  amplitudes  will  give  this 
phase  difference. 

With   Eqs.   (5) ,   the   equations  of  motion  for  the  blade 
become, 


Lift/span  =  -  a>2M  h   -  W2S  a   +  h  K  (IV-6a) 

b  o       b  o     Oh 


Moment/span  =  -  <o2S  h   -  w2I  a   +  a  K  (IV-6b) 

b  o       a  o     o  a 

The  natural  frequencies  of  twist  and  flexure  are  found  by 
assuming  that  each  mode  can  exist  separately,  in  the  absence 
of  exciting  forces.  Then, 

~7n  (IV-7a) 

h   b 
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<o       =    /K   /I 
a  a      a 


(IV-7b) 


K       =    Ma;2 
h  h 


(IV-8a) 


(IV~8b) 


and   therefore, 


icuT 
Lift/span   =    [-h    (o>2M      -    M  w2)    -    a  oj2S     ]e  (IV-9a) 

o  b  b    h  °         b 


Moment/span    =   [-hw2S      -    a    (w2I      -   I    w2)  ]e  (IV-9b) 

ob  o  a  a    a 


B.       TWO    DEGREE    OF    FREEDOM    FLUTTER    EQUATIONS 

The  method  used  to  investigate  flutter  in  this  report 
follows  the  procedure  and  notation  used  by  Garrick  and 
Rubinow  [Ref.  9],  who  derived  the  lift  force  and  moment 
equations, 

iwT      h 
L  =   -4/>b3u;2e        [    — °  (L      +    iL    )    +    a    (L      +    iL    )  1 

G/R  r  L    b         1  2  o      3  a     J 

(IV-10a) 

ia>T      h 
H  =   -4pb«a;2e        [  — 2  (M      +    iM    )     ♦    a     (M      +    iM    )  ] 

G/R  H  b         1  2  o  v    3  4' 

(IV-10b) 

where   L     was   defined   as   positive   down.   Equating  the 
G/R 

magnitudes  of   the   lift   from   Eqs.   (9a)   and   (10a)  ,   one 
obtains, 
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— 2  (L      +   iL      -  a   +  SI   X)     *   a    (L      +    iL      -    ux   )     =    0 
b        1  2  h  03  4         >a 

(17-11) 

and   by   equating   the   amplitude   of   the    moments, 

hQ 

—  (M      +    iM      -    /i.x  )    +    a    (M      +    1M      -    u.rz    +  &   X)     =    0 

b1  2a  03  n         r  a  h 

(IV-    12) 

where  the  non-dimensional  terms  are  defined  as, 

(L   *  iL  ) , (M   +  iM  )    -   real  and  imaginary  component 
12     12 

of  lift  and  moment  due  to  plunge. 

(L   +  iL  ) , (M   ♦  iM  )    -   real  and  imaginary  component 
3      4     3      4 

of  lift  and  moment  due  to  pitch. 

Mb 
u.         -   wing  density  parameter   (  — - —  ) 

Uphz 

£1         -   ratio  of  natural  frequencies  of  flexure 
h 

eu.  2 
and  torsion  p-(—~) 
a 

£1         -      the  quantity  /ir2. 

x    -   location  of  the  wing  center  of  gravity  measured 
a 

sb 
from  the  elastic  axis   ( — — ) 
\b' 

I 
r     -»   radius  of  gyration  of  the  wing   ( — — ) 

X    -   ratio  of  natural  frequency  in  pitch  to  the 

<on    2 
circular  frequency    (— -) 
ui 

Experimental   results   indicate   the    u        >    u>     in  modern 

a      h 

turbomachines  £Ref.  18]. 


h 
Eqs.   (11)   and   (12)  form  an  eigenvalue  problem  in  (— — ) 
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and  (a  )  which  has  a  non-trivial  solution  if  and  only  if  the 
o 

determinant   of   the   coefficient  matrix  vanishes.   Equating 

the  determinant  to  zero  yields  a   complex   function   in   the 

unknown    X,   each   component  of   which  must  be  identically 

equal  to  zero.  These  equations  are, 

H   a  X2    +   [ii    (L      -  fi)     +  SI    (M      -    ft  )  ]X    +    C      =0     (IV-13a) 
ha  a1  h3a  B 


(ft  L      ♦  ft   M    )  X    +    C      =0 

a  2  h    4  I 


(IV-13b) 


C      =   ia[  x    (M       +L)^(M      -ft)-Lr2-//.x2]+D 


1    a 


R 
(IV-13c) 


ci^[VV  V 


L    r2  J    *    D 
2    a  I 


(IV-13d) 


D      =LM      -    L   M      -    L    M      +LM 
R  13  3    1  2    4  4    2 


(IV-13e) 


D      =    L   U      -LM      +LM      -LM 
I  14  4    1  2    3  3    2 


(IV-13f) 


/X 


For  a  fixed  cascade 
geometry  and  interblade 
phase  angle,  Eqs.  (13) 
are  all  transcendental 
functions  of  the  reduced 
frequency  k,  or  (1/k) , 
more  accurately.  The 
solution  of  Eqs.  (13a) 
and  (13b)  is  found  using 
the  Theodorsen  method. 
/Re  {X}  and  /Im  {X}  are 
tabulated  and  plotted  versus  (1/k)  .  The  flutter  point  is 
defined  by  the  intersection  or  the  real  and  imaginary  curves 


/MxJ 


Figure    9.      Solution     plane     for 
the       Theodorsen       method. 


of  /X~ .  The  approximate  shape  of  these  curves  is  shown  in 
figure  9.  Since  (1/k)  appears  transcendentally  in  every 
term  of  Eqs.  (13),  the  curve  of  /Re {X}  is  not  exactly 
parabolic,  and  for  /Im [X}  ,  not  actually  linear. 


At    the   flutter  point,   the   non-dimensional   flutter 
frequency  is. 


%  =  /F 

and  the  non-dimensional  flutter  speed  is, 

UF   _  1  1 

cue     k/X 

k  is  defined  by  Eg.  (11-39)  . 


(IV-14a) 


(IV-14b) 


SINGLE  DEGREE  OF  FREEDOM  FLUTTER 


For  single  degree  of  freedom  oscillations  in  pitch,   the 
equation  of  motion  for  the  blade  is, 


Moment/span  =  I  a  ♦  aK 


(IV-  15) 


which  with  the  assumption  of  simple  harmonic  motion  becomes, 

itoT 


Moment/span    =   [-  cu2I    a      +    K   a    ]e 
a   o 


(IV-  16) 


The   mechanical   damping   of  the  system  is  introduced  by 

assuming  the  spring  opposing  the  displacement   in   pitch   to 

have   a   spring   constant  of  the  form   K  (1  +  ig  )  [Ref.  8  J. 

a       a 

The   natural    frequency   in      pitch      is      determined      as      before. 
Then,       equating      the      amplitude      of      this    expression    for    the 
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h 
moment  with  that  of  Eg.  (10a)  with  (— 2-)  =  0#  one  obtains, 

b 

a    (  X  -  1  )  +  M   =0  (IV-17a) 

a  3 


M   +  g  a   X  =  0  (IV-17b) 

U     a  a 


The  solution  method  is  the  same  as  for  the  two  degree  of 
freedom  problem,  and  the  non-dimensional  flutter  frequency 
and  speed  are  found  from  Egs.  (14).  A  similar  development 
is  required  to  obtain  the  equations  for  single  degree  of 
freedom  plunge.  However,  based  on  the  results  of  this,  and 
other  works,  plunge  oscillations  are  always  stable. 

For   both   single   degree   of   freedom  and   coupled 

th 
flexure-torsion  motion,  the  flutter  calculations  for  the  n 

blade   must   be  repeated  for  varying  interblade  phase  angles 

until  the  minimum  flutter  speed   is   found.   This   speed   is 

termed   the  critical  flutter  speed.  For  all  speeds  below  it, 

cascade  operation  will  be  stable. 

An  alternate  method  for  determining  the  onset  of  flutter 
is  to  compute  the  real  component  (the  in  phase  component)  of 
the  work  per  cycle  of  the  blade  on  the  freestream, 

Work/cycle  =  Re {  /  (M  ddt  +  L  hdt)}  (IV-  18) 

/    a        h 
Jo 

As   long   as  the  blade  does  positive  work  on  the  freestream, 

the  cascade  will  be   stable.   However,   if   the   aerodynamic 

damping   decreases  to  the  point  that  the  airstream  does  work 

on  the  blade,  flutter   will   occur.   For   single   degree  of 

freedom   pitch   oscillations,   with  zero  structural  damping, 

the  sign  of  the  imaginary  component  of  the  moment  determines 

whether   the   work/cycle   is   positive   or  negative,  and  can 

therefore   be   used   to   indicate   the   onset   of    flutter 
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[Ref.  4,18]. 


D.   DEFINITION  OF  TERMS 

The  complex  pressure  difference  at  each  grid  point  along 
a  blade  was  shown  by  Eg.  (111-66)  to  be  a  direct  function  of 
the  sonic  velocity  perturbation  amplitude.  Egs.  (111-76) 
then  gave  the  non-dimensional  lift  force  and  moment  per  unit 
span.  These  integrals  were  numerically  computed  using  a 
trapezoidal  integration  routine.  If  F  represents  the 
integrated  value  of  (m+1)  discret  values  of  f  (i) ,  then, 

F  =   f  (x)  dx  =  {-£f(1)  +  f(m+1)]  +  Vf(i)}Ax   (IV-  19) 


To  relate  the  non-dimensional  expressions  for  lift  force 
and  moment  as  derived  by  Garrick  and  Rubinow  to  those 
computed  using  the  Teipel  amplitude  functions,  the  following 
guantities  are  defined: 

hrt  iwT 

P(x,y,t)  =  [—  (PJ  +  a(P  )  ]e  (IV-20a) 

c    h      °  a 

Then  the  diraensionsal  lift  force  and  moment   per   unit   span 
are, 

1  1     h  ikt 

Lift  =  — />u2cL(x,y,t)  =  — pu2c[—  (L  )  +  a  (L  )  }  e 

(IV-20b) 


1  1      h0  ikt 
Moment  =  -pu262M  (x, y ,  t)  =  — pu2c2[—  (M  )  +  a  (M  )  ]e 

2  0  2   o   Le    h     o   a 

(IV-20c) 

L  ,  M  ,  L   and  M   are  all  complex  guantities.   The  subscript 
h    h    a       a 


indicates  the  perturbation  source.  Garrick  and  Rubinow 
defined  the  lift  force  as  positive  down.  Therefore,  equating 
Eqs.  (20b, c)  to  the  negative  of  the  corresponding  components 
of  Eqs.  (10) t    one  obtains, 


1   ho  ho 

—  [—  (L  )  ♦   (L  )  ]  =  [— -  (L   +  iL  )  +  a  L   +  iL  )  ] 
k*  c    h       a     Lb    1      2      03      4  J 


(IV-21a) 


2  K  K 

l-ZT   (M  )  ♦   (M  )  ]  =  [—  (K   +  iM  )  +  a  (M   *  iM  )  ] 

k2Lc    h        a      Lb    1      2      0   3      4  J 

(IV-21b) 

using   the  definition  of  reduced  frequency  from  Eq.  (11-39) . 
Hence, 


L   =  Re{L  } 

1    2k2     h 


(IV-22a) 


L   =  Im{L  } 

2    2k2      h 


(IV-22b) 


L   =  — Re{L  } 
3    k2     a 


(IV-22c) 


L   =  Im{L  } 

U        k2     a 


(IV-22d) 


M   =  Re  £M  } 

1    k2     h 


(IV-22e) 


2    k2      hJ 


(IV-22f) 


M   = Re{M  } 

3    k2      a 


(IV-22g) 
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2 

M   = Im{M  }  (IV-22h) 

The   computer   program,  which  is  listed  after  the  appendices 

and  is  explained  in  Section  V,  computes  the   non-dimensional 

lift  force  and  moment  defined  by  Eqs.  (20b, c)  for  each  blade 

of  the  cascade.  These  values  are  then  modified  by  Eqs.   (22) 

th 
so   that  the  flutter  analysis  can  be  made  for  the  n    blade. 
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COMPUTER  PROGRAM 


A.   INTRODUCTION 

The  finite  difference  equations  derived  in  Section  III-B 

and   the  flutter  equations  of  Section  IV  were  programmed  in 

FORTRAN  IV  and  run  on  an  IBM  360/67  computer.  The   computer 

program   listing  for  the  two  degree  of  freedom  case  of  pitch 
and  plunge  is  given  following  Appendix  B. 

The  computational  procedure  is  the  same  for  pitch  or 
plunge,  the  only  difference  between  the  two  being  the  form 
of  the  flow  tangency  equation.  The  two  degree  of  freedom 
solution  is  simply  a  superposition  of  the  individual  cases 
of  pitch  and  plunge. 

It  should  be  noted  that  the  subscripts  of  the  four  grid 
points  of  the  computational  molecule  shown  in  figure  3  are 
distinct  from  the  notation  used  for  subscripting  the 
perturbation  quantities  in  the  computer  program.  There,  the 
subscripts  indicate  whether  or  not  a  perturbation  quantity 
at  a  grid  point  is  above  or  below  a  surface,  and  if  the 
perturbation  is  due  to  oscillation  in  pitch  or  plunge. 

In  the  computer  program,  all  perturbation  quantities 
have  an  alpha-numeric  subscript  (two  numbers  and  a  letter) 
and  an  index.  The  index  identifies  the  grid  point.  The 
letter  in  the  subscript  tells  whether  the  component  is  real 
(R)  or  imaginary  (I) ,  while  the  numbers  indicate: 

22  (pitch)  or  44  (plunge)  perturbation 

quantity  is  at  a  grid  point  which  is  not 
on  the  lower  surface  of  a  blade  or  wake 
slip  plane. 
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33  (pitch)  or  55  (plunge)  perturbation 

guantity  is  at  a  grid  point  located  on 
the  lower  surface  of  a  blade  or  slip 
plane. 

Thus,  V    is  the  complex  component  of  the   perturbation   in 
551 

normal   velocity,   due  to  plunge,  on  the  lower  surface  of  a 

blade  or  slip  plane.  The  subscript  of  the  perturbation 
guantity  solved  for  in  this  section  merely  indicates  the 
real  or  imaginary  component  of  that  guantity  at  grid  point 
22. 


PflOGBAM  CONSTRAINTS 


Figure  10.   The  compatible  stagger  angle. 

The  computer  program  is  quite  general  in  its  ability  to 
solve  a  given  cascade  flow  problem.  The  one  geometric 
constraint  that  must  be  satisfied  by  the  user  is  Eg.  (II-1) , 
subsonic  leading  edge  locus  of  the  cascade  blading.  The 
computer  program  requires  that   the   leading   edge   of   each 
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blade  be  exactly  on  a-  grid  point.  Therefore,  the  input 
stagger  angle  is  modified  so  that  STGR,  shown  in  figure  10, 
is  a  non-zero  integer  multiple  of  DSTSTR  aft  of  the  initial 
Mach  wave.  Then  the  'compatible'  stagger  angle  is  determined 
and  given  in  the  output.  The  amount  by  which  this  procedure 
alters  the  input  cascade  geometry  is  a  function  of  the 
gridf ineness,  Mach  number,  solidity,  and  distance  between 
blades.  See  also  V-G. 

The  only  other  restrictions  are  on  vector  size.  Internal 
checks  are  made  at  the  beginning  of  a  run  to  insure  that  the 
program  dimensions  are  large  enough  for  the  input  geometry 
and  Mach  number.  If  the  storage  area  is  exceeded,  program 
termination  occurs.  Of  course,  the  program  is  only  valid  for 
M  >  1,  and  in  the  range  where  linearized  theory  is  valid. 

The  vector  dimension  requirements  for  a  given  cascade 
and  Mach  number  can  be  determined  prior  to  running  the 
program.  The  vectors  used  for  storage  of  all  perturbation 
quantities  and  grid  point  locations  along  a  right-running 
Mach  line  must  be  dimensioned  an  integer  value  greater  than, 

1  1 

-€  tan  (5  )  {(N   -1)[tan(/3)  -  cot  (  a  )]  +  -}+  1  (V-  1) 

2  bid  d 

where, 

€  •  •  •  •  grid  fineness  ratio 

N     •  •  *  •  number  of  blades 
bid 

d     o   •  •  •  distance  between  blades 

and  a.    and  f3   have  been  defined  previously. 

The  program  marches  down  successive  right-running  Mach 
lines  from  one  grid  point  to  the  next.  The  location  of  a 
grid   point,   along   the   current   Mach   line,   at  which  the 
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perturbation  quantities  have  just  been  computed  is  IHAVEP. 
The  next  grid  point  is  labeled  IHAVEP+1.  In  the 
subroutines,  these  same  two  grid  points  are  labeled  J  and 
I  respectively.  Then  IHAVEP+1  =1  is  the  index  of  the 
storage  location  within  a  vector  array  of  a  perturbation 
quantity  or  grid  point  location  at   F     of   figure   3.   The 

left-running   Mach  line  which  passes  through  grid  points  F 

1 1 

and  F    is  a  line  of  constant   IHAVEP   =   J,   and   the   next 
21 

left-running  Mach  line  is  a  line  of  constant  IHAVEP+1  =  I. 


The  last  grid  point  in  the  flow  field  is  the  first  wake 
grid  point  aft  of  the  last  blade  in  the  cascade.  Along  the 
left-running  Mach  line  which  passes  through  this  grid  point, 
IHAVEP  =  LIMIT.  When  the  computation  reaches  this  Mach  line, 
the  program  switch  variables  are  reset,  and  the 
computational  procedure  jumps  to  the  next  right-running  Mach 
line.  Thus,  perturbation  quantity  indices  are  from  1  to 
LIMIT+1,  while  IHAVEP  goes  from  0  to  LIMIT  along  a  Mach 
line. 

The   perturbation   quantity  vectors  are  single  dimension 

arrays  in  order  to  reduce  required  core  space.   However,   as 

the   computational  molecule  requires  perturbation  quantities 

from   two   different   Mach   lines,   it    is   necessary   to 

temporarily   switch   and   store   the   perturbation   quantity 

values  as  scalars  upon  entrance  to   each   subroutine   before 

solving   for   'new'   values   at   F   .    Figure   11  shows  the 

22 

storage  locations  of  a  representative  perturbation   quantity 

vector   as  it  is  in  the  program  on  entrance  to  and  exit  from 
a  subroutine. 


ENTRANCE 


IHAVEP-I 


IHAVEP+I=I 


Perturbation    quantities 
from   the    current    right- 
running     Mach    line. 


Perturbation     quantities 
from  the    previous    right- 
running     Mach    line. 


EXIT 


IHAVEP=J 


IHAVEP+I=I 


Figure  11.   Typical  perturbation  guantity 
vector  storage  on  entrance  to  and  exit 
from  a  subroutine. 


The  computational  procedure  in  subroutines  GENFPT, 
NEWBLD,  TOP,  BOTTOM  and  HAKE  is  the  same.  This  procedure 
is: 

Step   1.   The   incomming   perturbation  quantities  at 

grid  point  IHAVEP=J  on  the  current  Mach  line  (F    of   figure 

21 

3)   are   stored  as  the  scalar  quantities  having  suffixes  NEW 

(U2RNEW,  V4INEW,  etc) .  Stored  in  the  vector  on  the  left  at 
location  IHAVEP  are  the  perturbation  quantities  from  the 
previous  Mach  line  (F   )  . 

Step   2.   The  constants  K12R,  K12I,  etc.,  defined  in 
Section   III-B,   are   computed   using    the   perturbation 
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quantities   at   grid   points   F    and  F   .  Now,  perturbation 
u  11       21 

quantities  from  the  previous  Mach  line,  stored  at   location 
IHAVEP  (vector  on  the  left)  ,  are  no  longer  required. 

Step  3.  Perturbation  quantities  from  the  previous 
grid  point  of  the  current  Mach  line  (U2RNEW,  V4INEW,  etc.) 
are  stored  in  the  vector  at  position  IHAVEP=J  (vector  on 
right) ,  which  removes  from  memory  the  perturbation 
quantities  from  the  previous  Mach  line  at  grid  point  IHAVEP. 

Step  4.  The  perturbation  quantities  at  IHAVEP+1  =  I 
are  computed  and  stored  as  the  scalar  quantities  U2RNEW, 
V4INEW,  etc.,  and  the  subroutine  is  exited. 

This  procedure  is  slightly  different  in  subroutines   BOTTOM, 

WAKE   and   NEWBLD   because  perturbation  quantities  for  lower 

surface  points  are  computed.  When  a  lower  surface  point  has 

been   computed   by   one   of   these   subroutines,  the  integer 

variable  ICO  is  set  equal  to  one.  This   flagged   value   then 

causes   GENFPT,   which   will   always   be  the  next  subroutine 

called,  to  use  lower  surface  perturbation  quantities  in   the 

computational  molecule  at  grid  point  F    for  the  calculation 

21 

of  perturbation  quantities  at  grid  point  F 

22 

The  last  grid  point  on  all  right-running  Mach  lines  will 
lie  in  either  the  general  flow  field  or  the  wake  of  one  of 
the  blades.  When  this  grid  point  is  reached,  after  returning 
from  the  appropriate  subroutine  to  MAIN,  the  perturbation 
quantities  just  computed  are  stored.  Then  the  program  jumps 
to  the  next  right-running  Mach  line.  The  program  terminates 
flow  field  calculations  on  the  first  grid  point  in  the  wake 
aft  of  the  last  blade. 

Subroutine  COEF  stores  the  pressure  difference  across  a 
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blade  resulting  from  pitch  or  plunge.   These  values  are  then 

used  to  compute  the  lift  and  moment  (about  the  elastic  axis, 

x  )  for  each  blade.  In  order  to  be  able   to   identify   these 
o 

pressure    differences   with   the   appropriate   blade,   the 

information  is  stored  in  a  two-dimensional  matrix.  The  first 
dimension  identifies  the  blade,  while  the  second  dimension 
specifies  the  maximum  number  of  grid  points  along  the  blade. 
In  order  to  minimize  the  required  core  size,  it  was 
desirable  not  to  have  the  first  dimension  equal  the  total 
number  of  blades  in  the  cascade.  Instead,  a  smaller  matrix 
was  used  cyclically  for  all  of  the  blades.  However,  this 
required  that  a  restriction  be  placed  on  the  number  of 
blades  that  any  right-running  Mach  line  can  cross. 

After  the  first  wake  point  aft  of  a  blade  is  reached, 
C0EF2  (alternate  entry  point  for  subroutine  COEF)  is  called, 
and   the   pressure   difference  at  the  end  point  of  the  blade 

(x    =  1)  linearly   extrapolated.   Then   the   dimensionless 
rel 

lift   force   and   moment,   as   given   by   Egs.  (111-76) ,  are 

computed  using  a  trapezoidal   integration   technique.    Once 

the   lift   and  moment  have  been  computed  that  portion  of  the 

matrix  storing  information  relative  to   that   blade   may   be 

reused  for  another  blade.   The  maximum  number  of  blades  that 

a  right-running  Mach  line  may   cross   is   ICROSS,   which   is 

specified  in  subroutine  INPUT  by  a  data  statement.   In  order 

not  to  exceed  the  storage  space  allocated  for  these   vectors 

storing  pressure  and  moment  distributions,  the  right-running 

Mach  line  which  passes  through  the  leading  edge  of  the  (n   + 

th 
ICROSS   +   1)     blade   must  pass  aft  of  the  first  wake  grid 

th 
point  of  the  n    blade.    Also,   there   must   be   less   than 

MAXPTS  grid  points  along  a  blade,  including  the  extrapolated 
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point  at  x    =1.   ICEOSS  and  MAXPTS  are  the  dimensions   of 
•         rel 

the  vectors  storing  the  pressure  and  moment  distributions. 


Letting  J  |  denote  the  chopped  integer  value  of  a  real 
number,  the  restriction  on  blades  crossed  is  satisfied  if 

1 

ICROSS{|d[tan  {Q  ) -cot  (  a  )  ]/  Ax|     +€}    >    I +    2    1 

Ax 

(V-  2) 

The  right  hand  side  of  Eg.  (2)  gives  the  total  number  of 
grid  points  along  the  blade  (which  must  be  less  than 
MAXPTS) .  x  is  the  step  size  along  a  streamline  (DSTSTR) , 
and  is  equal  to 

Ax  =  2d[cot(a)  ]/€  (V-  3) 

All  terms  were  previously  defined  in  Eg.  (1) . 

One  final  dimension  that  must  not  be  exceeded  is  the 
total  number  of  blades,  which  can  not  be  greater  than 
MAXBLD.  MAXPTS  and  MAXBLD  are  specified  by  a  data  statement 
in  INPUT  in  the  same  manner  as  ICROSS.  Integrated  lift 
forces  and  mcments  are  stored  by  blade  number.  From  figure 
11,  it  can  be  seen  that  apparent  blades  are  'seen*  by  the 
computer.  Depending  on  the  input  parameters  (low  Mach 
number  and  high  solidity),  the  computer  may  'see'  several 
apparent  blades.  Therefore,  MAXBLD  should  be  at  least  two  or 
three  greater  than  the  number  of  blades  in  the  cascade.  The 
required  core  size  is  insensitive  to  the  magnitude  of 
MAXBLD. 

Thus  ,  the  constraints  and  restrictions  on  program  use 
are  given  by  Eg.  (II-1)  and  Eqs.  (V-1)  -  (V-3)  and  the 
comment  on  the  maximum  number  of  blades  above.  All  of  the 
constraints  are  checked  internally  by  the  pr  >gram  in  INPUT, 
and  an  appropriate  error  termination  message  given  if  one  or 


more  of  the  constraints  is  not  met. 


PROGRAM  FLOW  AND  CONTROL 


With  the  exception  of  the  main  program  and  subroutine 
COEF,  all  subroutines  operate  in  a  straightforward  manner  to 
solve  the  finite  difference  equations  developed  in  Section 
III-B.  Hence,  they  will  be  only  briefly  discussed. 

MAIN  controls  the  program  flow  as  consecutive 
right-running  Mach  line  are  traversed.  Integer  switch 
variables,  listed  in  Appendix  B,  are  used  in  the  the  five 
primary  logic  tests.  These  are,  in  the  order  encountered, 

(1)  If  (  LCOUNT  =  0  )  grid  point  is  on  the  initial 

right-running  Mach  line 

(2)  If  (  IHAVEP  =  LCOUNT  )  grid  point  is  on  a  blade 

or  in  a  wake 

(3)  If  (IHAVEP  =  0  )  grid  point  is  on  the  initial 

left-running  Mach  line 

(4)  If  (  IHAVEP  =  NUMPI*NSTPTS  )  and  (JLINE  =  MAXI  } 

grid  point  is  on  the  leading  edge  of  a 
new  blade 

If  (2)  was  satisfied,  then, 


(5)   If  (  IHAVEP  >  LIMW  )  grid  point  is  in  the  wake 

To  clarify  program  operation  and  the  use  of  these  switch 
variables,   the  sample  case  shown  in  figure  12  is  discussed. 

For  M  =  1.3   and   the   geometry   given,   the   quantities 
listed  in  figure  12  are  computed  prior  to  beginning  any  flow 
field  calculations.    •   is   the   compatible   tagger   a] 
computed   by   the   program.  Subroutine  INITAL  calculates  the 
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perturbation  quantities  at  the  leading   edge   of   the   first 

blade   and  initializes  IHAVEP  =  1,  JLINE  =  0  and  LCOUNT  =  0. 

MAIN  then  computes  NEWDST  and  LIMIT  as  shown  in  the   figure, 

and  initializes 

LIMW  =  NEWDST  =  10 

MAXI  =  NGRDFN  +  NSTPTS  =  10 

NUrl  =  0 

NUMOLD  =  0 

NUMPI  =  1 

OLDL  =  0 

ICO  =  0 

Now  the  loop  which   calculates   the   entire   flow   field   is 

entered. 


IHAVEP  =  0 

— (D 

Assumed 

/        blade 

INPUT 
M      =     I.3 
€      =     6 
/?    =     62.8 
d      =    .33824 
Nb.d   =    2 
COMPUTED 
Ax   =  .09365 
Or      =   I.356 
&'      =      62. 71 
NSTPTS    =    4 
NEWDST    =    10 
LIMIT      =     15 


IHAVEP=  LIMIT 


Figure    12.      Sample   Cascade   Problem 


IHAVEP      =      0    means   that    the    grid    point    is  on   the    initial 

left-running    Mach    line,    while      IHAVEP       *      0       ,.       LCOUNT      =      0 
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indicates  that  the  grid  point  is  on  the  initial 
right-running  Mach  line.  In  either  case,  subroutine  MACHLN 
is  called  to  compute  the  the  perturbation  quantities 
associated  with  each  grid  point.  For  computations  everywhere 
except  along  the  initialright-running  Mach  line,  the 
temporary  shuffling  of  perturbation  quantities  must  be  done 
prior  to  computing  the  values  at  the  new  grid  point. 

After  computing  down  the  initial  Mach  line  to   IHAVEP   = 

LIMIT,   the   program   arrives  at  the  first  grid  point  on  the 

second  right-running  Mach  line  with  IHAVEP  =  0  and  LCOUNT   = 

JLINE   =   1.  All  other  integer  switch  variables  are  the  same 

as  before.  At  the  second  grid  point  on  this  Mach  line  (point 

1   in   figure   11),   IHAVEP  =  LCOUNT  =  1,  which  results  in  a 

call  to  TOP/BOTTOM  (since  IHAVEP  <  LIMW) .   Prior  to   calling 

TOP,   ICO  is  set  equal  to   1  so  that  the  next  time  GENFPT  is 

called  (at  2)  the  perturbation   quantities   from   the   lower 

surface   of   the   blade   will   be   used  in  the  computational 

molecule  for  grid  point  F   .  After  return  from  BOTTOM,  C0EF1 
21 

(alternate   entry  to  subroutine  COEF)  is  called  to  store  the 

pressure  difference  just  computed  across  the  first  blade  at 

this  grid  point.  Then  IHAVEP  is  incremented  by   1. 

For  the  remaining  grid  points  along  this  Mach  line, 
since  all  of  the  logical  tests  will  be  false,  GENFPT  will  be 
called  by  default.  Again,  ICO  =  1  indicates  that  a  lower 
surface  point  is  to  be  used  for  perturbation  quantities  at 
grid  point  F    in   the   computational   molecule.   Prior   to 

exiting   GENFPT,   ICO   is   reset   equal   to   zero.   For   the 

remainder  of  the  grid  points  on  the  Mach  line.  GENFPT  will 
be  called  until  IHAVEP  =  LIMIT.  The  perturbation  quantities 
computed  at  each  grid  point  are  stored  in  vector  locations 
IHAVEP+1.  When  IHAVEP  =  LIMIT,  the  integer  switches  are 
reset  to 

IHAVEP  =  0 
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LCOUNT  =  LCOUNT+1  =  2 

JLINE  =  JLINE* 1  =  2 

ICO  =  0 

ICO   is   reset   to   zero  here  to  allow  for  the  case  when  the 

Mach  line  terminates  on  a  grid  point  in  the  wake  (as  at   3)  , 

since  a  call  to  WAKE  causes  ICO  =  1. 

The  computational  procedure  remains  the  same  through  the 

tenth   Mach  line  (along  A)  .  MACHLN  is  called  whenever  IHAVEP 

=  0,  and  TOP/BOTTOM  are  called  when  IHAVEP  =  LCOUNT.   GENFPT 

is   called  otherwise.  At  the  completion  of  computation  along 

(A)  ,  the  integer  switch  variables  are  reset  to 

IHAVEP  =  0 
LCOUNT  =  10 
JLINE  =  10 
ICO  =  0 

and  the  other  integer  variables  are  as  before, 

LIMW  =  10 

NSTPTS  =  a 
NUMPI  =  1 
MAXI  =  10 

The  computational  procedure  down  Mach  line  (B)  is   the   same 

until  the  grid  point  at  (4) ,  where, 

IHAVEP  =  (NUMPI)  (NSTPTS)  =  4 

and 

JLINE  =  MAXI  =  10 

Satisfaction   of   this  logical  test  envokes  a  call  to  NEWBLD 

to  compute  the  perturbation  quantities  at  the  leading   edge 

of  the  second  blade.  NEWBLD  also  also  increments  NUM  so  that 

NUM  =  1  on  return  .  Next,   C0EF1   is   called   to   store   the 

pressure   difference   at   the  first  grid  point  on  the  second 

blade.  Next, the  integer  switches  are  reset  to, 

ICO  =  1 

LCOUNT  =  IHAVEP  +  NGRDFN  =  10 

JLINE  =  IHAVEP  =  4 

IHAVEP  =  IHAVEP  +1=5 

OLDL  =  IHAVEP  =  5 

NUHPI  =  NUM  +1=2 

MAXI  =  NGRDFN  +  (NUMPI) (NSTPTS)  =  14 

NUM  =  NUM  -1=0 

Computation   down  (B)  continues  until  the  grid  point  at  (5)  , 

with  GENFPT  being  called  for  all   intervening   grid   points. 

For   the   first  grid  point  after  (4),  ICO  =  1,  and  was  reset 


equal  to  zero  in  GENFPT.   At  (5)  ,  IHAVE  =  LCOUNI  =   10,   and 

TOP/BOTTOM   are   called.   COEF1   is   called   on   return  from 

BOTTOM,  and  GENFPT  is  called  for  the  remaining   grid   points 

until  IHAVEP  =  LIMIT.   Then  the  integer  switch  variables  are 

reset  to 

IHAVEP  =  0 
NUM  =  NUMOLD  =  1 
ICO  =  0 
JLINE  =  JLINE  +1=5 

Since  more  than  one  blade  has  now  been  encountered,  NUMOLD  > 

0,  and  it   is   necessary   to   reset   the   following   integer 

switches   so   that   the   grid   point   on   blade   two   can  be 

identified. 

LCOUNT  =  OLDL  =  5 

LIMW  =  (NUMOLD)  (NSTPTS)  +  HEWDST  =  14 

OLDL  =  OLDL  +1=6 

The   computation   now    starts   over   at   the   initial 

left-running  Mach  line  (Mach  line  C)  ,  where  MACHLN  is  called 

for   IHAVEP   =   0.    Thereafter,  until  the  grid  point  at  (6) 

where  IHAVEP  =  LCOUNT  =  5,  GENFPT  is  called.   At  (6),  ICO  is 

set   egual   to  one  and  TOP/BOTTOM  are  called.  Next,  C0EF1  is 

called  to  store  the  pressure  difference  just  computed.   Upon 

return   from  COEF1,   since   NUM   >   0,   the  integer  switch 

variables  must  be  reset  so   the   grid   point  on   the   upper 

surface   of   the   next   blade  (or  wake,  in  this  case)  can  be 

identified.  Thus, 

NUM  =  NUM  -1=0 

LCOUNT  +  NGRDFN  =  11 

LIMW  =  LIMW  -  NSTPTS  =  10 

GENFPT  is  then  called  to  compute  the  perturbation  quantities 

at  grid  points  until  the  program  reaches  (7) ,  where  IHAVEP  = 

LCOUNT   =   11.   This   time,  IHAVEP  >  LIMW  ,  which  means  that 

WAKE  should  te  called  instead  of  TOP/BOTTOM.  Again,   ICO   is 

set   equal   to  one   before  calling  WAKE.   After  return  from 

WAKE,  COEF2  is  called.   As  this  is  the  first   time   that   it 

has  been  called  for  this  blade  (i.e.,  this  is  the  first  grid 

point  in  the  wake  aft  of  blade  one),  the  integration  of   the 

pressure   distribution   for   the   lift   and   moment  is  done. 


Since  this  is  the  first  call  to   C0EF2   aft   of   blade   one, 

IWRITE   =.  1   on  return  to  the  main  program  and  the  pressure 

distribution  for  this  blade  can  be  output,  if  desired.  Next, 

a   check   is   made   to  determine  whether  or  not  this  was  the 

last  grid  point  in  the  flow  field,  or  if  IHAVEP  =  LIMIT.   If 

not,  the  calculation  loop  is  reentered,  with  ICO  =  1.  GENFPT 

is  called  to  compute  the  remaining  grid  points  until   IHAVEP 

=  LIMIT,  where  the  integer  switches  are  reset  to 

IHAVEP  =  0 
NUM  =  NUMOLD  =  1 
ICO  =  0 
JLINE  =  JLINE  +1=6 

and  as  NUMOLD  >  0, 

LCOUNT  =  OLDL  =  6 

LIMW  =  (NUMOLD)  (NSTPTS)  +  NEWDST  =  14 

OLDL  =  OLDL  +1=7 

The  computational  procedure  down  Mach  line  (D)  is  the 
same  as  for  (C) .  The  call  to  COEF2  at  (8)  will  result  in  an 
immediate  return  to  the  main  program  with  IWRITE  =  0  as  this 
will  be  the  second  time  (more  accurately,  not  the  first 
time)  that  COEF2  has  been  called  for  blade  one.  The 
computation  proceeds  in  the  same  manner,  down  successive 
right-running  Mach  line,  until  reaching  point  (9)  on  Mach 
line  (E) ,  where, 

IHAVEP  =  (NUMPI)  (NSTPTS)  =  8 
and 

JLINE  =  MAXI  =  14 
which  initiates  a  cal,l  to   NEHBLD   for   the   assumed   blade. 
Termination  of  the  flow  field  calculations  is  indicated  when 

NUM  >  NUMBLD  -1=1 
and 

IHAVEP  =  LIMIT 
at  point  (10).  As  (10)  is  also  the  first  wake  grid  point  aft 
of  blade  2,  COEF2  is  called  to  compute  the  lift   and   moment 
for   the   second   blade   before   terminating   the  flow  field 
calculations. 
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D.   SUBROUTINE  COEF 

Subroutine  COEF  is  a  multiple  entry  subroutine  used  to 
store  the  pressure  difference  distributions  for  each  blade 
as  it  is  encountered,  to  extrapolate  for  the  pressure 
difference  at  the  trailing  edge  of  each  blade,  and  to  use 
these  distributions  to  compute  the  lift  and  moment  due  to 
pitch  and  plunge.  If  the  pressure  distribution  for  a 
particular  blade  is  required,  output  is  obtained  by  issuing 
a  write  statement  after  return  from  C0EF2  in  the  main 
program. 

COEF  is  the  first  section  of  this  subroutine.  Called 
after  return  from  INITAL,  it  initializes  all  arrays  used  in 
the  subroutine  and  stores  the  pressure  difference  for  the 
first  grid  point  on  the  first  blade.  The  pressure 
difference  is  then  used  to  compute  the  incremental  moment 
about  the  elastic  axis. 

C0EF1  is  called  after  every  call  to  T0P/3OTT0M  or 
NEWBLD,  and  stores  the  pressure  difference  with  the  real 
component  referenced  to  the  blade  initial  position,  and  the 
imaginary  component  referenced  to  the  mean  position,  blade 
pitching  leading  edge  up.  Refer  to  Eqs.  (111-70) 
(111-75) .  The  pressure  differences  are  then  used  to  compute 
the  incremental  moment  resulting  from  each  discrete  pressure 
difference. 

C0EF2   is   called   after   every  call  to  WAKE.   The  first 

time  this  section  is  called,  it  initiates  the  calculation  of 

th 
the   lift   force   and   moment   for   the   n    blade,   after 

extrapolating  for  the  pressure  difference   at   the   trailing 
edge   of   the   blade.    At  the  same  time,  the  integer  switch 
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ID{n)  is  set  equal  to  one.  Subsequent  calls  to  C0EF2  for  the 
n     blade   wake   result   in  an  immediate  return  to  the  main 


program,  after  setting  IWRITE  =   0.   IWRITE   is   an   integer 

switch    variable    used    to   indicate   when   the   pressure 

distribution  for  the  entire  blade  is  complete.  For   a   given 

blade,   IWRITE  =  1  only  after  the  first  call  to  C0EF2,  where 

the  pressure  difference  at  the  trailing  edge  of  the  blade  is 

determined.   At  all  other  times,  IWRITE  =  0.   After  the  lift 

force  and  moment  have  been   computed,   the   section   of   the 

th 
matrix    used   for   storage  of   the   n     blade   pressure 

distribution  can  be  reused  for  another  blade.    Upon   return 

th 
to   MAIN,  the  pressure  distribution  for  the  n   blade  can  be 

output,  if  desired,  if  IWRITE  =  1. 


E.   OTHER  SUBROUTINES 

The  remaining  subroutines  used  to  calculate  perturbation 

quantities   in   the   flow   field   are   TOP,   BOTTOM,  GENFPT, 

NEWBLD,  tfAKE,  and  MACHLN,  all  called  by   the   main   program. 

Each   subroutine   solves   the   applicable   finite  difference 

equations  derived  in  section  III-B.  Subroutine  WAKE  uses  IBM 

subroutine  SIMQ,  a  simultaneous  equation  solving  routine,  to 

compute  the  perturbation  quantities  at   F   .    The   solution 

22 

procedure   in   each   subroutine   is   the   four   step  process 

described  earlier  in  this  section,  and  it  is  carried  out 
first  for  pitch,  and  then  for  plunge  in  each  subroutine. 
Each  subroutine  defines  .several  constants,  a  tew  of  which 
are  common  to  both  solutions. 


Subroutine  INPUT  reads  the  input  data,  computes  the 
compatible  stagger  angle,  and  the  geometric  cascade 
parameters  used  by  other  subroutines.  Here  the  input  data  is 
checked   to  insure  that  available  core  size  is  not  exceeded. 

Subroutine  INITAL  uses  the  equations   derived   for   the 
initial  Mach  lines  to  compute  the  perturbation  quantities  at 
the  leading   edge   of   the   first 
blade.  In  addition,  JLINE,  LCOUNT 
and  IHAVEP  are  initialized. 


DSTSTR 


Figure    13.       Characteristic 
mesh     geometry. 


Finally,  subroutine  COMPXY 
computes  the  grid  point 
coordinates  in  the  characteristic 
mesh  in  terms  of  the  the  geometry 
computed   in   INPUT   and  shown  in 

figure  13.  The   calculation   procedure   for   points   on   the 
right-running  Mach  line  is 

X  (I)  =  X  (I)  +  HDSTRL 

Y(I)  =  Y(I)  -  TRNGLH 

The  first  grid  point  on  a  right-running  Mach  line  is 
X(1)  =  X  (1)  +  HDSTRL 

Y(1)  =  Y(1)  +  TRNGLH 

Actually,  the  only  coordinate  used  by  the  program  is  X (I) . 


Subroutine  FLUTER  is  called  for  two  degree  of  freedom 
flutter  solutions  when  the  input  variable  IFLUT  =  1,  and 
uses  the  equations  derived  in  Section  IV,  where  the  single 
degree  of  freedom  equations  were  given  also.  Flutter 
investigation  is  quite  expensive  in  terms  of  computer  time, 
especially  if  the  program  is  designed  to  iterate  to  an 
actual  flutter  speed,  subject  to  some  error  tolerance.  To 
minimize  the  'first  blade  effect*,  the  flutter  calculatio] 
must  be  conducted  for  a  relatively  high  numbered  blade  where 


the  magnitude  of  the  lift  forces  and  moments  have  more 
nearly  approached  an  assymptotic  value.  This  must  be  done 
because  an  actual  turbomachine  does  not  have  a  'first1  blade 
and  the  aerodynamic  forces  on  the  blades  are  assumed  to  be 
periodic.  To  determine  a  critical  flutter  speed  for  a  given 
cascade  geometry  and  Mach  number,  the  minimum  flutter  speed 
must  be  determined,  which  requires  an  iteration  on  the 
interblade  phase  angle  as  well.  Since  each  iteration 
through  the  flutter  routine  requires  a  recomputation  of  the 
entire  flow  field,  it  is  apparent  that  computational  time 
mounts  rapidly. 

In  an  attempt  to  reduce  the  amount  of  computer  time 
required  for  the  flutter  investigation,  the  computer  program 
is  set  up  to  compute  and  output  the  roots  of  the  flutter 
determinant  versus  the  reciprocal  of  the  reduced  frequency 
(1/k) .  These  results  are  then  plotted  by  hand,  using  the 
Theodorsen  method  to  find  a  flutter  solution. 

c 1/k)  is  incremented  by  an  amount  STEP,  an  input 
variable.  If,  for  a  given  value  of  (1/k) ,  the  real  or 
imaginary  component  of  the  solution  is  negative  the  program 
will  terminate,  as  the  dimensionless  flutter  frequency  would 
be  imaginary.  However,  if  the  discriminant  of  the  quadratic 
for  Re  {X}  is  neqative,  then  (1/k)  is  first  decreased  by 
STEP,  and  then  increased  by  one-fifth  that  amount.  The  flow 
field  is  then  recomputed,  and  FLUTER  called  again.  Note  that 
this  procedure  assumes  that  the  reduced  frequency  used  for 
the  previous  solution  produced  a  positive  real  root.  As  the 
real  component  of  the  flutter  determinant  is  quadratic,  a 
negative  discriminant  indicates  that  the  value  of  (1/k)  is 
past  the  vertex  of  the  parabolic  curve  of  /Re  {X}  in  the 
solution  plane  (see  figure  9,  Section  IV)  .  Each  time 
through  the  flutter  subroutine,  the  square  root  of  the  three 
roots  of  the  flutter  determinant  and  the  corr  ponding  valu< 
of  (1/k)  are  output.   Because  of  the   manner   in   which   the 
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case  of  a  negative  real  root  is  handled,  and  also  because 
the  minimum  flutter  speed  occurs  at  the  smallest  value  of 
(1/k) ,  the  program  should  initially  be  given  a  reduced 
frequency  on  the  order  of  three  or  greater.  In  all  cases 
investigated  in  this  report,  all  three  of  the  flutter 
determinant  roots  were  positive  for  reduced  frequencies  of 
that  magnitude. 

F.   INPUT  AND  OUTPUT 

Subroutine  INPUT  reads  all  input  data  from  two  Namelist 
input  records,  NAM1  and  FLTR.  The  input  variables  for  NAM1 
are, 

NGRDFN  -  the  grid  fineness  ratio  (e)  .  The  number 
of  increments  along  a  Mach  line  between 
two  blades. 

FSTRMN  -  the  freestream  Mach  number. 

REBFRQ  -  the  reduced  frequency,  defined  in  Eg. 
(11-39) . 

XSUBO   -  the  elastic  axis  position  (x  ) 

o 

TNWDST  -  vertical  distance  between  adjacent 
blades. 

STGANG  -  the  input  stagger  angle  (/5)  .  The 
compatible  stagger  is  computed  in   INPUT. 

FAZE   -  the  interblade  phase  angle  (in  degrees) . 

NUMBLD  -  the   number   of   blades   (N    )   in    the 

bid 
cascade. 

IPRES1  -  number  of  the  blade  user  desires  specific 

information     about,     i.e.,  pressure 

distribution,   lift   force   and  moment 
values,  etc. 

IPRES2  -  same  as  IPRES1. 


Variables   in  Namelist  FLTR  are  all  used  in  the  flutter 
routine.  They  are, 

RMUU    -  the  blade  density  parameter  (fx)  . 

XSUBA   -  location   of   the   center    of   gravity 
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measured   from  the  elastic  axis,  positive 
aft  (xa). 

RSUBA   -  the  blade  radius  of  gyration  (r  ) . 

WHWA   -  ratio  of  natural  frequency  in   plunge   to 
the  natural  frequency  in  pitch  (— ) 

STEP    -  increment  for  O/k)   . 

IFLUT   -  set   equal   to  one  (IFLUT=1)  if  a  flutter 
solution  is  desired. 


The  input  data  for  the  sample  program  discussed  in  part 
C  of  this  section  would  be: 

&NA£J1 

NGRDFN  =  6, 
NUMBED  =  2, 
FAZE  =  50.  , 
XSUBO  =  0.5, 
TNWDST  =  .33824, 
STGANG  =  62.8, 
IPRES1  =  1, 
IPRES2  =  2, 
SEND 
&FLTR 

RMUU  =  500., 
XSUBA  =  0.. 
RSUBA  =  0.5, 
HHWA  =  0.70  7, 
STEP  =  0.05, 
IFLUT  =  0, 
&END 

Each   card  begins   in   column   2,   and  the   comma   follows 

immediately  after  the  last  digit.  Spaces  between   the   last 

digit   and  the   comma   are   interpreted  by  the  computer  as 
zeros. 

Output  can  be  modified  at  the  users  discretion.  With  the 
program  as  listed  after  Appendix  B,  pressure  distributions, 
lift  forces  and  moment  and  the  magnitude  and  phase  of  the 
complex  lift  and  moment  for  each  blade  can  be  output.  In 
addition,  if  appropriate  'write1  statements  are  put  in  each 
of  the  subroutines,  the  entire  flow  field  may  be  printed. 
This  option  produces  sizeable  quantities  of  output,  and  is 
useful  primarily  for  program   debugging.    The   output 
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flutter  solution  runs  was  described  in  Section  V-D. 


G.   GENERAL 

Because  the  program  computes  a  compatible  stagger  angle, 
the  cascade  geometry  actually  used  in  the  computer  program 
may  be  slightly  different  than  that  input  by  the  user.  For 
large  values  of  the  grid  fineness  ratio,  e  ,  the  differences 
are  insignificant.  However,  it  is  normally  preferable  to  use 
a  coarser  grid  size  to  minimize  computer  time.  As  a  result, 
if  a  specific  cascade  is  to  be  input,  then  an  iteration  on  e 
should  be  done  to  find  the  smallest  value  of  e  which 
produces  the  compatible  stagger  angle  that  is  in  closest 
agreement  with  the  input  stagger  angle.  If  this  procedure  is 
not  followed,  the  cascade  computed  for  will  be  slightly 
different  than  desired,  and  therefore,  the  reflection  points 
on  the  blades  will  be  different  than  expected. 

The  computer  program  requires  approximately  175K  bytes 
of  storage  (G  compiler)  or  135K  bytes  of  storage  (H 
compiler)  when  run  on  an  IBM  360/67  computer.  Compile  time 
is  approximately  twenty-five  seconds.  Actual  CPU  time  per 
run  is  dependent  on  the  input  geometry,  grid  fineness  and 
number  of  blades.  Sample  times  are  given  for  each  of  the 
results  discussed  in  Section  VI. 

Initial  program  checkout  and  verification  can  be 
accomplished  in  several  ways.  By  setting  the  reduced 
frequency  equal  to  zero,  the  well  known  result  for  the  lift 
on  a  flat  plate  at  unit  angle  of  attack  is  obtained  (for  the 
first  blade) .  Also,  tabulated  results  appearing  in  the 
report  by  Garrick  and  Rubinow  [Ref.  9]  for  single  blades 
require  less  than  five  seconds  of  computer  time  to  generate 
using  a  grid  iineness  of  100  or  less. 
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VI.  RESULTS 


Initial  work  was  done  for  single  degree  of  freedom 
torsion.  All  cases  discussed  in  this  section  were  run  with 
zero  structural  damping  (g   =   0)  .    Pressure   distributions 

and  flutter  solutions  were  computed  and  compared  with  Snyder 

[Ref.  19],  who  used  a  finite  difference  procedure  developed 
by  Verdon  [Ref.  23].  Brix  [Ref.  2]  showed  that  considerable 
discrepancy  existed  between  the  pressure  distributions  given 
by  Snyder  and  those  computed  by  the  method  of 
characteristics,  especially  in  the  region  of  reflection 
points.  As  stated  by  Verdon  [Ref.  25f  p24], 

The  finite  difference  marching  procedure  used  .  .  .  does 
not  recognize  pressure  discontinuities  as  such,  but  only 
approximates  them  by  sharp  changes  in  the  pressure 
distribution.  This  numerical  approximation  also  gives 
rise  to   the   ripples   which   appear   in   the   predicted 

Sressure     distributions      downstream     of      the 
isco n tin ui ties. 

The  flutter  comparisons  shown  in  figures  13  and  14  are 
for  cascades  C,  D  and  E,  described  below. 

Cascade   C      Cascade   D      Cascade   E 
B      =    67         $      =    60         B     =   70 
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The  entry  in  the  Time  column  is  the  actual  CPU  time,  in 
seconds,  required  for  one  iteration  through  the  flow  field 
for  a  ten  blade  cascade.   Figures  begin  on  page  91. 

All   three  cascades      have   a   solidity   of   1.33   and  a 


mid-chord  pitch  axis  (x  )  .   The  grid  fineness,  e  ,  was  varied 
o 

in  order  to  be  able  to  match  the  geometry  of  each  cascade  as 

closely  as  possible  (refer  to   V-G)  ,   while   maintaining   at 

least   twenty-five   discrete   grid   points  along  each  blade. 

Both   results  assumed   zero   structural   damping   (g  ) .   In 

addition,   the   method   of   characteristics  solution  used  an 

assumed  value  of  500  for  the  blade  density  parameter  //.,  and 
a   radius   of  gyration,  r  ,  of  0.5.  The  value  for  ft  was  felt 

to  be  reasonable  for  modern  turbomachinery  blading  [Ref.  5]. 

These   figures,   14   and  15  show  the  reduced  freguency  which 

occurs   at   the   minimum   flutter   speed,    k  ,and   the 

crit 

interblade   phase   angle  at   that  point,  8     ,  versus  Mach 

crit 

number.   Method  of  characteristics  results  are  based  on   the 

tenth  blade.  Agreement  is  not  exact.  However,  in  light  of 
the  previous  comment  on  the  computational  scheme  used  by 
Snyder,  the  agreement  is  felt  to  be  reasonable. 


In  figures  16  to   20,   similar   results   are   shown   for 

Cascade  C  with  elastic  axis  as  a  parameter.  Again,   /x  =   500, 

and  r   =  0.5. 
a 

Figures   21   to   24  present  flutter  stability  boundaries 

for  Cascade  C.  These  are   plots   of   elastic  axis   position 

versus   Mach   number  for  which  M   =0.  The  shaded  portion  of 

4 

each  graph  is  unstable,  i.e.,  M   <  0.  The  analysis  was   made 

4 

for  30°  increments  in  the  interblade  phase  angle,  8 .  Figures 

21  and  22  showed  instabilities   at   only   the   values   of 
shown.   In   figure   23,  the  range  of  interblade  phase  angles 
over  which  instabilities  existed  increased,   and   in   ligure 
24,   instabilities   occurred   for   all  values  of    examined. 


Thus,  these  four  figures  show  that  as  the  reduced  frequency 
is  increased  (1/k.  decreased)  ,  cascade  operation  becomes  more 
stable. 

When  conducting  flutter  investigations,  it  is  important 
to  find  the  highest  value  of  reduced  frequency  for  which  the 
cascade  becomes  unstable,  as  this  will  produce  the  lowest 
flutter  speed  for  a  given  geometry  and  interblade  phase 
angle.    Figure   25  shows  the  variation  of  M   with  (1/k) .  As 

can  be  seen,  if  the  step  size  in  (1/k)   is   too   large,   the 

first   point   at  which  H   <  0  could  be  missed  (note  inset  in 

the  figure).   As  a  result,  the  reduced  frequency  at  which  H 

4 

=   0   at   B   would  result  in  an  incorrect  (too  high)  flutter 

speed.  Also,  figure  25  shows  the   condition   where   M    just 

4 

crosses   the  axis.   When   the   phase   angle   was   slightly 

decreased,  the  curve  of  M   was  not  quite  tangent  to  the  axis 
4 

at   A,   but   crossed   the   axis   near  B.  Figure  26  shows  the 

resulting  jump  in  the  flutter  speed  and  reduced  frequency  at 

flutter. 


A  limited   amount   of   investigation   was   conducted   to 

determine   the   dependence  of  flutter  speed  on  blade  number. 

Although  M   is  quite  oscillatory  for  the   first   few   blades 
4 

and   does   not   approach  an  asymptotic  value  until  after  the 

tenth  or  fifteenth  blade,  preliminary  results,  based  on  a 
limited  number  of  cases,  indicated  that  the  flutter  speed 
and  frequency  were  fairly  constant.  An  example  is  shown  in 
figure   27  for  Cascade  D  with  M  =  1.7,  8  =  -45  ,  and  e   =  25. 

Comparison  with   results   by   Verdon   [Ref.  25]   for   an 
infinite   cascade   was   made  using  the  two  degree  of  freedom 


computer  program  listed  at  the   end   of   this   report.    The 

cascades   A   and   B  referred  to  in  the  figures  are  shown  and 

described  in  figure  28.   Pitch  results  are   always   for   the 

fifteenth   blade,  and  plunge  results,  for  the  fourteenth.  In 

figures  29  through  49,  the  time  dependence  is   in   terms   of 

the   compressible   reduced   frequency   as  defined  by  Verdon, 

M 
where   k   =   k     .    Corresponding   to   Verdon's   notation, 
c     M2-1 

Im {M  }   >   0   implies  unstable  operation  for  the  pitch  case. 

Plunge   oscillation   was   always   found   to   be   stable.   In 

addition,   figure   29   shows  convergence  of  the  norm  of  M 

(norm   =   |M  J   =  /H  ®M~  )   versus   blade   number   for   four 
a  a  a 

different   values   of  8  .   These  phase  angles  are  all  in  the 

region  for  which  Verdon  was  unable  to  obtain  a  solution. 
The  curves  are  shown  as  being  continuous  for  clarity  of 
presentation,  but  results  were  only  obtained  at  discrete 
integer  values  of  blade  number,  of  course.  Each  curve 
required  approximately  nineteen  minutes  of  computer  time. 

Figures  30  through  33  are  plots  of  the  real  and 
imaginary  component  of  the  pressure  distribution  acting  on 
the  upper  and  lower  surfaces  of  the  blades  for  the  two 
cascades.  Figures  34  through  37  show  the  pressure  difference 
across  the  blade  (lower  -  upper)  for  two  different 
interblade  phase  angles.  Agreement  was  excellent. 

Figure  38  defines  the  symbols  used  in  figures  39  through 
49,  where  lift  and  moment  loops  are  obtained  by  plotting  the 
real  and  imaginary  components  of  the  lift  force  and  the 
moment  with  the  interblade  phase  angle  as  a  parameter.  The 
method  of  characteristics  was  able  to  calculate  solutions 
throughout  the  range  of  interblade  phase  angles,   0   <   8  < 


360.   Again,  excellent  agreement  was  obtained. 

In  figures  30  through  49,  results  for  cascade  A  required 
48-50  seconds  of  CPU  time  per  fifteenth  blade  solution,  and 
for  cascade  B,  100-110  seconds  per  fourteenth  blade 
solution. 

The  last  figure,  50,  shows  a  flutter  solution  for 
cascade   A  assuming  the  following  parameters:  jj.    -    500,  r   = 

.707. 
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Figure  14.   Reduced  frequency  at  minimum  flutter  speed  versus 
Mach  number  '  a =500, r  =0.5  ). 
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Figure  15.   Interblade  phase  angle  at  minimum  flutter  speed 
versus  Mach  number  (/x=500,r  =0.5  )T 
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Figure  17.  Reduced  frequency  (kcrit)  an<3  interblade  phase 
angle  (Scrit)  at  minimum  flutter  speed  versus  Mach  number 
(/x  =  500,ra=0.5  ). 
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Figure  18.  Reduced  frequency  (kcrit)  and  interblade  phase 
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Figure  19.  Reduced  frequency  (kcr^t)  and  interblade  phase 
angle  (8cr-[t)  at  minimum  flutter  speed  versus  Mach  number 
(  ^=500,ra=0.5  ). 
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Figure  20.  Reduced  frequency  (kcrlt)  and  interblade  phase 
angle  (Scr-Lt)  at  minimum  flutter  speed  versus  Inch  number 
(/i=500,ra=0.5  ). 
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Figure  21.   Flutter  boundaries  versus  elastic  axis  position 
(shaded  portion  is  unstable) . 
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Figure  22.   Flutter  boundaries  versus  elastic  axis  position 
(shaded  portion  is  unstable). 
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Figure  23.   Flutter  boundaries  versus  elastic  axis  position 
(shaded  portion  is  unstable) . 
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Figure  24.   Flutter  boundaries  versus  elastic  axis  position 
(shaded  portion  is  unstable) . 
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Figure  25.   Oscillation  of  the  imaginary  component  of  the 
moment  due  to  pitch  (M^)  with  reduced  frequency. 
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Figure  26.   Variation  of  flutter  speed  and  reduced  frequency 
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Figure  27.   Non-dimensional  flutter  speed  and  reduced  frequency 
at  flutter  versus  blade  index  ( yx=500,r  =0.5  ). 
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Figure  28.   Cascade  A  and   B  geometry. 
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Figure  29.  Convergence  of  the  norm  of  the  moment  due  to  pitch 
oscillation  with  blade  index  as  a  function  of  the  interblade 
phase  angle. 
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Figure  30.  Comparison  of  the  pressure  distribution  resulting  from 
pitch  oscillation  of  the  fifteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   A) . 
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Figure  31.  Comparison  of  the  pressure  distribution  resulting  from 
plunge  oscillation  of  the  fourteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   A) . 
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Figure  32.  Comparison  of  the  pressure  distribution  resulting  from 
pitch  oscillation  of  the  fifteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade  B) . 
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Figure  33.  Comparison  of  the  pressure  distribution  resulting  from 
plunge  oscillation  of  the  fourteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   B) . 
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Figure  34.  Comparison  of  the  pressure  difference  distribution  (lower  - 
upper)  resulting  from  pitch  oscillation  for  the  fifteenth  blade  with 
the  infinite  cascade  results  of  Verdon  (Cascade   A). 
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Figure  35.   Comparison  of  the  pressure  difference  distribution  (lower 
upper)   resulting  from  plunge  oscillation  of   the  fourteenth   blade 
with  the  infinite  cascade  results  of  Verdon  (Cascade  A). 
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Figure  36.  Comparison  of  the  pressure  difference  distribution  (lower  - 
upper)  resulting  from  pitch  oscillation  for  the  fifteenth  blade  with 
the  infinite  cascade  results  of  Verdon  (Cascade   B) . 
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Figure  37.  Comparison  of  the  pressure  difference  distribution  (lower  - 
upper)  resulting  from  plunge  oscillation  of  the  fourteenth  blade 
with  the  infinite  cascade  results  of  Verdon  (Cascade   B) . 
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Figure  38.   Definition  of  symbols  used  in  figures  39  through  49. 
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Figure  39.  Comparison  of  the  moment  coefficients  due.  to  pitch 
oscillation  for  the  fifteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade  A) . 
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Figure  40.  Comparison  of  the  moment  coefficients  due  to  pitch 
oscillation  for  the  fifteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   A) . 
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Figure  41.  Comparison  of  the  moment  coefficients  due  to  pitch 
oscillation  for  the  fifteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   B) . 
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Figure  42.  Comparison  of  the  moment  coefficient  due  to  pitch 
oscillation  for  the  fifteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   B) . 
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Figure  43.  Comparison  of  the  lift  coefficients  due  to  plunge 
oscillation  for  the  fourteenth  blade  with  the  Infinite  cascade 
results  of  Verdon  (Cascade  A) . 
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Figure  44.  Comparison  of  the  lift  coefficients  due  to  plunge 
oscillation  for  the  fourteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade  .  A) . 
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Figure  45.  Comparison  of  the  lift  coefficients  due  to  plunge 
oscillation  for  the  fourteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade  A) . 
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Figure  46.  Comparison  of  the  lift  coefficients  due  to  plunge 
oscillation  for  the  fourteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   B) . 
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Figure  47.  Comparison  of  the  lift  coefficients  due  to  plunge 
oscillation  for  the  fourteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   B) . 
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Figure  48.  Comparison  of  the  lift  coefficients  due  to  plunge 
oscillation  for  the  fourteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   R) . 
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Figure  49.  Comparison  of  the  lift  coefficients  due  to  plunge 
oscillation  for  the  fourteenth  blade  with  the  infinite  cascade 
results  of  Verdon  (Cascade   B) . 
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Figure  50.  Non-dimensional  flutter  speed  and  reduced  frequency  at 
flutter  versus  interblade  phase  angle  for  Cascade  A  (/x=500,ra= 
0.5  ). 
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VII.  CONCLUSIONS 


A  method  of  characteristics  procedure  has  been  developed 
for  the  solution  of  the  complete  linearized  unsteady 
aerodynamics  of  a  flat  plate  cascade  having  a  subsonic 
leading  edge  locus  in  a  supersonic  flow  field.  Pressure 
distributions  were  computed  and  the  resulting  lift  forces 
and  moments  used  to  investigate  flutter  boundaries  and 
stability  limits  for  a  two-dimensional  finite  cascade,  which 
can  have  arbitrary  stagger  angle  and  interblade  phase  angle. 

Agreement  with  the  infinite  cascade  results  by  Verdon 
were  excellent.  However,  further  verification  and  comparison 
with  new  methods  being  developed  by  others  [Refs.  12,20,26] 
should  be  conducted.  Of  primary  importance  would  be  the 
development  of  a  periodic  solution  applicable  for  the 
infinite  cascade. 
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VII. 


CYLINDRICAL  SHELL 


This  developement  is  closely  based  on  a  previous  papar 
by  Platzer,  Brix  and  Webster  [Ref.  17]. 

A.   PROBLEM  EORfiULATION 

Consider  the  supersonic  flow  of  an  inviscid,  adiabatic 
gas  through  a  circular  cylinder  whose  axis  is  aligned  with 
the  freestream,  and  whose  surface  is  rigid  upstream  of  the 
origin.  The  cylinder  geometry  is  shown  in  figure  51. 
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m  =  I 
n  =  2 


Figure  51.   Cylinder  Geometry 

Downstream   of   the  origin,  a  flexible  panel  of  length  L 

is   assumed   to   undergo   small  amplitude   vibrations.   The 

governing   equation   for   this  problem   is   the  linearize!, 

unsteady   velocity    potential  equation    in   cylindrical 
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coordinates,    $(X,R,9,T), 

11  M  1 

(1-H2)<j,         +    <j,         +    -4,       +    _  $         -    2-<> <(>         =    0 

XX  RR         R    R  R2    99  C    XT         c2    TT 

(VIII-    1) 
The   linearized   form    of    the    pressure   coefficient   is, 

2  2 

C      = <{> 4>  (VIII-    2) 

p  u2    T        uo    X 

The  linearized  boundary  condition  en  the  cylinder  inner 
surface  is 

dh  dh 

4>    |  =   — -  +    u  —  ,  0<X<L  (VIII-   3) 

R|r=R0         dT  o^X  ~ 

The   boundary   condition  at  r  =  0    is  discussed  separately. 

In  Eg.  (3),  h  =  h(X,9,T)   is   the   equation   for   the   panel 

surface.  One  other  condition  that  must  be  satisfied  is  the 
Sommerfeld  radiation  condition;  i.e.,  waves  must  propagate 
away  from  the  source  of  disturbance.  For  a  derivation  of  the 
above  three  equations,  refer  to  Bisplinghof f ,  Ashley,  and 
Halfman,  "Aeroelasticity"  [Ref.  1]. 

To  non-dimensionalize  the  above  equations,  lengths  are 
referred  to  flexible  cylinder  length  and  time  to  flexible 
cylinder  length  divided  by  the  freestream  velocity.   Then, 

1  1 

(1-M2)4>         +    <|>         +   -<|>       +   — 1>         -    2M2<{>         -    M24>         =    0 
xx  rr        r   r        r2    99  xt  tt 

(VIII-    4) 


C      =    -2[<J)      +    <f>    ]  (VIII-    5) 

p  t  x 


subject    to, 
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ah     ah 

t>     I      =  —  +  —  ,     0  <  x  <  1  (VIII-  6) 

r|r=R    at   ax        * 

The  non-dimensional  shell  deflection  is  assumed  to  be  of  the 
form, 

ikt 
h(x,9,t)  =  Z(x)cos(n0)e  (VIII-  7) 


where, 


k  =  —  (VIII-  8) 


and   n    is   the   circumferential   mode   number.    Then  the 
perturbation  potential  must  be  of  the  form, 

ikt 
<Mx,r,6,t)  =  0(x,r)cos(n0)  e  (VIII-  9) 

and  the  resulting  equations  are, 

1  n* 

(1-112)0     +  0  +  -0         +   (  k2M2 >  )0    - 

xx    rr   r  r  r2 

i2kH0   =  0  (VIII-10a) 


C   =  -  2(>   +  ikjzf]  (VIII-10b) 

P        x 

C   is  the   pressure   coefficient   amplitude.    The   boundary 

P 
conditions  are, 


0    I      =  Z   +  ikZ  ,  0  <  x  <  1            (VIII-11a) 
r  |r=Ro     x 

The   boundary  condition   at  the   axis   is  dependent  on  the 

circumferential  mode  number  n.  Thus, 

0    \           =  0  ,    for   n  even                 (VIII-11b) 
r  |r=0 


C  I     =  0  ,    for   n   odd  (VIII-11c) 

p|r=0 
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Egs.   (10)   and   (11)   are  the  basic  equations  which  will  be 
used  in  the  cylindrical  shell  analysis. 

B.   METHOD  OF  CHARACTERISTICS 

For  M  >  1,  Eg.  (10)  is  hyperbolic  and  has  real 
characteristics  which  satisfy  the  ordinary  differential 
equation 


(M2-1) d*r  -  d*x  =  0 


(VIII-  12) 


For  ds  the  differential  arc  length  along   a   characteristic, 
d*s  =  M2d*r,  and, 


dx 
ds 

= 

/M2^1 

M 

dr 
ds 

= 

1 

±  — 
H 

(VIII-13a) 


(VIII-13b) 


The  derivative  along  a  characteristic   of   an   arbitrary 
function  F(x,r)  is, 


dF    dFdx   dFdr 
ds   dxds   drds 


(VIII-    14) 


Then, 


dF         /M2-1  1 

F      =   —    =      F      +  —  F 

i         ds,  tt         x         M    r 


(VIII-15a) 


dF        /M2-1  1 

2         ds,  M        x        Mr 


(VIII-15b) 


where      ds        and      ds        are   the   differential    arc   lengths   along 

1  2 

the  left  and  right  characteristics  respectively.   The  second 
(cross)  derivatives  are, 
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-F    =  -[  (1-M2)F    *  F   ]  =  -  F 
*2   M        xx    rr       2 


(VIII-  16) 


From  Eg.  (15) , 


(F   +  F  J 


x    2/M2-1   i     2 


(VIII-17a) 


F   =  —  (F   -  F  ) 
r    2   »     2 


(VIII-17b) 


Now,   with   Eqs.   (16)   and   (17) ,   the   basic  equations  are 
written, 


:  T~(  0     -   J*   )  ♦  (  *2  '  -TT7  )* 


12   '2i    2rM    »     2 

kM 


r2M2 


:(    0*0) 
/M2-1      1       2 


(VIII-18a) 


t     =  -  (0     -   0  )    =   Z   +  ikZ 
r    2   i     2     x 


(VIII-18b) 


C   =  -  2[0      +  ikjzf] 
P         x 


(VIII-18C) 


C.   FINITE  DIFFERENCE  EQUATIONS 


1.   General  Field  Point 


The   finite   difference   equation  for  the  points  P  and  S 
along  the  right-running  Mach  line  is, 


<VS)  -*a<P»  "  ££<*,  "  "«»  +AS(  ^    'S^)0 


kMAs 
-  i7r==r(  rf   +  J*   ) 


(VIII-  19) 
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Note   that  in  this  and  subsequent  equations,  ds   =  ds   =  As. 

12 

For   points   Q   and   S    along   the   left-running    Mach   line. 


As  n2 

(S)     -   ft    (Q)     =  — -  (rf      -   0   )  +  As  (   R2    -  )  0 

*  2rM      i  2  M2r2 


kMAs 
i-^==r(    0+t) 


/M2-1 
The  computational   molecule    is   shown   in   figure  52. 


(VIII-  20) 


Figure  52.   Characteristic   grid  with 
computational  molecule. 


The  values  for  0,  0  and  0  are  averages  between 
adjacent  grid  points  along  the  respective  Mach  lines,  so  for 
the  left-running  Mach  line. 


1 

=  -C0x(Q)    +  0X  (S)  ] 


(VIII-21a) 


1 


-102(Q)     +    02{S)  ] 


(VIII-21b) 
and   0 (S)     is    obtained    by    integrating    along    the   Mach    line. 
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0(S)     =  0(Q)     +  -C^(Q)  ♦  tl  (S)  ]AS 


(VIII^21c) 


0   =  "If  (S)    +  J*(Q)  ] 


Similarily,  along  the  right-running  Mach  line, 


(VIII-21d) 


0l    =   ~10 t(P)  +  0t(S)   3 


(VIII-22a) 


02  =  -C^2(P)  +  02(S)  ] 


(VIII-22b) 


0<S)  =  >*(P)  +-[^2(P)  +^z(S)]As 


(VIII-22c) 


^  =-[p(P)  +  0(S)  ] 


(VIII-22d) 


Egs.   (21)  and  (22)  are  substituted  into  Egs.  (19)  and  (20)  . 
The  result  is  a  system  of  two  eguations  for  &    (S)  and  0    (S) 


of  the  form, 

K0    (S)  +  B0    (S)  =  C 


(VIII-23a) 


D0    (S)  +  E0    (S)  =  F  (VIII-23b) 

1  2 

Eg.  (23a)  applies  to  the  right-running  Mach  line,  where 

(VIII-24a) 


As       kMAs 
A  =  1 +  i- 


UMr(S)     2/M2-1 


As 
B  =  -  (  k2 

4Mr  (S) 


n2     As  2      kHAs 

)  ( )    +  i   ; 

M2r2(S)     2       2/P^l 


(VIII-2Ub) 


135 


C    =  0,(2)1 


As 
4Mr  (P) 


+     (    k2    - 


n2  As    2  kMAs 

)  ( — )       -    i  ] 

M2r2  (S)  2  2/M2^! 


As  kMAs 


4Mr  (P)  2/M2-1 

n2  n2  As 

♦   0(P)[  (   k2 -  )    +    (   k2 )  ](-— ) 

M2r2(P)  M2r2(S)  2 

(VIII-24c) 
Eg.     (23bJ    applies    to    the  left-running    Mach    line,    with, 


As 
4Mr  (S) 


E    =   -    1 


+     (    k2    - 


As 


4Mr(S)  2/M2-1 


n2     As  2 

kMAs 

M2r2  (S)     2 

2t/M2_rT 

(VIII-25a) 

kMAs 

L .   ... 

(VIII-25b) 

F    =   tf    (Q)[ 


As 


4Mr  (Q) 


-    (   k2 


n2  As   2  kMAs 

)  ( — )      +   i     .     .     ] 

M2r2(S)  2  2/PrT  J 


As  kMAs 


4Mr(Q)  2/M2-1' 


n2  n2  As 

*  (Q)  t  (  *2  ^  „a  9//xl  )  +  (  *2  -  Mg  rrr,   )  3  (-7) 


M2r2  (Q) 


M2r2(S)  2 

(VIII-25c) 


Egs.     (23)    are   used    to   compute  0  and      its      derivatives      at      a 
general   flow   field    point. 

2«       Initial   Mach    Line 

Along      the   initial   right-running   Mach   line,    0   =  0     = 
0,    and   Eg.     (23a)    reduces   to, 


(S)     =   C/A 


(VIII-26a) 
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and  C  simplifies  to. 

As       kMAs 

C    =  0    (P)[1  +  -  i— ;    .1  (VIII-26b) 

*         4Mr  (P)     2//M2:rT 


3.   Inner  Surface  Of  Outer  Cy_linde£ 

The   flow   tangency   condition  prescribes  the  normal 
velocity  at  r  =  R  .  Thus, 


M 
0    j      =  -O  (S)  -  0    (S)  ]  =  Z   +  ikZ        (VIII-27a) 
r|r=R0    2   * v  '     2        x 

For  the  initial  point  on  the  panel,  0     =   0,  and 


2  I 

0  (Sj  =-£Z   +  ikZ]  (VIII-27b) 

1       M   x        |x=0 

Thereafter,  the  equations  for  the   left-running   Mach   lines 
are  used  to  solve  for  0,    0      and  0    .  Hence, 


2 
(S)  =  [F  +  — (Z   +  ikZ)E]/(D  +  E)  (VIII-28a) 

M   x 


2 
(S)  =  [F  -  -(Z   +  ikZ)D]/(D  +  E)  (VIII-28b) 

M   X 


C   =  -  210      +  ik£]  (VIII-28C) 

P         x 

Then,  0      and  0      are  used  to  compute  0      and  C  . 
12  x       p 

4.   Boundary,  Condition  at  the  Axis 


The   boundary  condition  at  r  =  0  was  approximated  r>y 
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assuming  a  very  thin  solid  cylinder  with  radius  r  =  R   along 

i 

the   axis,  extending  the  length  of  the  cylinder.  The  form  of 

the  boundary  condition  applied  at  r  =  R   is  dependent  on  the 

i 

circumferential  mode  number   n.  For  even  values  of   n, 

0  =  0     n  even  (VIII-29a) 


Thus, 


(S)  =  C/(A  +  B)  =  0     (S) 


(VIII-29b) 


For   odd   values  of   n,  the  prescribed  boundary  condition  is 

C   =0.  For  the  steady  case  (k  =  0) ,  this  condition   reduces 

P 
to  0        =   0.  Then  applying  the  boundary  condition  with  Eqs. 
x 

(23a)  and  (17), 

C  I      0  =  0    (S)  +  iM(S)     n  odd  (VIII-30a) 

p  r=R,       x 


1 
0(S)     =  *(P)  ♦~J>1(P)  +^2(S)]AS 


(VIII-30b) 


0    (S)  = 


2/M2-1   i 


10    (S)  +  0    (S)  ]  =  -  ik^(S)     (VIII-30C) 


A^(S)  +  B02(S) 
Defining, 
8 


M        As 

.       ■■■■  +  ik  (— ) 
2v^2=T        2 


(VIII-30d) 


(VIII-31a) 


X  =  ik[0(P)  +  (— )02(P)  ] 


(VIII-31b) 
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Then, 


y  MB 

J*     (S)     =    [     C    +    B -£-]/[     A    -    -7===-    ] 
*  S  2/M2-1S 


(VIII-31c) 


0     (S)     =    [     C    -    A0     (S)y  ]/B 


(VIII-31d) 


Now      0(S)  ,    $    (S)  ,    and   0    (S)     are   computed    using    Egs.     (17)    and 
x  r 

(22a)  . 

The   initial   right-running  Mach  line  is  assumed  not 
to  reflect  from  the  assumed  inner  cylinder. 

D.   GENERALIZED  FORCES  AND  THE  PRESSURE  COEFFICIENT 

Having  obtained  the  pressure  distribution  over  the  inner 

surface  of  the  outer  cylinder  using  the  pressure  coefficient 

amplitude,  the  generalized  aerodynamic  force  Q     is   found 

mr 

from, 


2   ■  -  /C  (x)yp   (x)dx 

mr    2  /  pm    r 


where, 


Z(x)  =  Ea  <//  (x) 
3      3    3 


(VIII-33a) 


(VIII-33b) 


and   C   (x)   is  the  pressure  distribution  resulting  from  the 
pm 

th  th 

m   axial   mode   deflection   and  ^  (x)   is   the   r     axial 


deflection   mode.   Thus,   Q    would  be  the  generalized  force 
x  2 

due  to  an  axial  deflection  of  at   least   two  modes   of   the 
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form, 

Z(x)  =  A  ^  (x)  +  A  ^  (x)  (VIII-34a) 

11  22 

and, 

Q    =  —  /c   (x)V'  (x)  dx  (VIII-34b) 

12    2/  Pi     2 

The  above  equations  were  programmed  in  Fortran  IV  and 
run  on  an  IBM  360/67  computer.  For  200  grid  points  along  the 
flexible  cylinder  length,  total  run  time  was  approximately 
twelve  to  fifteen  seconds  when  using  a  precompiled  source 
deck.  The  program  is  listed  after  the  Two  Degree  of  Freedom 
cascade  listing. 


Pressure  coefficients  and  generalized  aerodynamic  forces 
were  computed  for  sinusoidal  deflection  modes, 

Z(x)r  E^  (x)  =  sin(JTrx)  ,    j  =  1,2,  •  •  •   (VIII-  35) 
3 


These  cases  were  previously  studied  by  Widnall  and 
Dowell,  [Ref.  27],  who  derived  a  solution  of  the  full 
linearized  equation  using  a  Laplace-transform  technique. 

For  figures  53  through  55  the  grid  fineness  (  number  of 
increments  along  a  Mach  line  between  the  inner  and  outer 
cylinder  )  was  adjusted  as  a  function  of  the  difference  in 
radii  (or  R/L)  to  maintain  200  points  along  the  unit 
cylinder  length.  The  inner  radius  was  kept  fixed  at  0.0002 
while  the  outer  radius  was  varied  from  0.1  to  5.0. 
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For  low  values  of  the  circumferential  mode   number    n  , 

pressure     distributions    were    excellent,    with    sharp 

discontinuities  evident  at  reflection  points  on   the   outer 

cylinder.    However,   for    n   >   2,  these  reflection  points 

exhibited  a  marked  oscillation  in  pressure   whose   amplitude 

and   duration   were   dependent   on  the  inner  radius  and  grid 

fineness   used.    By    arbitrarily    varying    these    input 

parameters,   the   oscillations  could   be  minimized,  but  not 

entirely  removed.  Over  a  small  range  of  inner  radii,  0.01  to 

0.0002,    overall   effect   on   Q     was   relative  small  with 

mr 

differences    of   less   than   fifteen   percent   in    most   cases.      For 

R/L      large      enough      that      no      reflection      occurred,       Q        and 

mr 

pressure   distributions   were    independent   of    inner   radius  size 

and  grid  fineness  used.  Further  investigation  is  being 
conducted  to  develop  stability  criteria  as  a  function  of 
Mach   number,    grid   fineness,    and   R/L. 


Figure    53.      Generalized   aerodynamic    force   Q;qr   versus    radius/length 
ratio   as   a   function  of    the   circumferential  mode  number    (Ri=.0002). 
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M  =  l.l 


Widnall    and    Dowell  (Ref.27) 

°    Present     study 

R/L  =   l.O 
n  =    o 


M  =  l.2 

■ 

M=  JZ 

i          i          i          i 

.0        .1        .2        .3 


.4       .5        .6        .7       .8        .9        1.0 
K 


Figure  54.   Generalized  aerodynamic  force  Qui  versus  reduced 
frequency  as  a  function  of  Mach  number  (Ri=.0002). 
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Figure  55.   Generalized  aerodynamic  force  Q12R  versus  radius/length 
ratio  as  a  function  of  the  circumferential  mode  number  (Ri=.0002). 
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F.   CONCLUSIONS 

The  method  of  characteristics  solution  of  the  full 
linearized  unsteady  velocity  potential  equation  has  been 
developed  to  compute  the  pressure  distributions  and 
generalized  forces  for  harmonically  oscillating  cylindrical 
shells  with  arbitrary  radius-to-length  ratios,  axial  and 
circumferential  mode  numbers,  supersonic  Mach  numbers  and 
reduced  freguency.  In  addition,  the  method  is  easily 
adaptable  to  include  arbitrary  (non-sinusoidal)  boundary 
conditions.  Moreover,  it  is  felt  to  be  more  easily  applied 
than  the  complex  integration  techniques  required  for  the 
Laplace-transform  solution. 


APPENDIX.  A 

TWO  DEGREE  OF  FREEDOM  CASCADE 

The   following   variables   appear  in  the  logic  statements  of 
the  program. 

ICO  used  in  subroutine  GENFPT.  When  equal  to  one, 

it  indicates  that  the  previous  grid  point  was 
a  lower  surface  (blade  or  wake)  point.  ICO  is 
always  set  equal  to  zero  before  exiting  from 
GENFPT. 

ID  integer  array  with  each  element  corresponding 

to  a  cascade  blade.  Normally  zero,  elements 
are  set  equal  to  one  the  first  time  C0EF2  is 
called  for  each  blade. 

IFLUT  an   input   variable.   Set   equal   to   one   if 

flutter  solutions  are  not  desired;  zero 
otherwise. 

IHALT  a   termination  flag.  When  equal  to  one,  X,  as 

determined  from  the  real  or  imaginary 
component  of  the  flutter  determinant,  is 
negative. 

IHAVEP  the  number  of  the  computation  step  along  a 
right-running  Mach  line.  Set  equal  to  zero  at 
the  initial  left-running  Mach  line,  it  is 
incremented  until  it  is  equal  to  LIMIT. 

IPRES1,IPRES2  input  variables  identifying  blades  for  which 
specific  information  is  desired 

ITER  iteration   counter   for   the   number  of  times 

through  the  flutter  routine. 

ITO  an   information   flag.   If  other  than  zero  on 

return  from  FLUTER,  the  discriminant  of  the 
quadratic  for  the  real  roots  of  X  is  less 
than  zero. 

ITflMAX         maximum  number  of  flutter  solution  iterations 
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IWEITE 


LIMIT 


NGBDFN 


NEWDST 


allowed  before  program  termination. 

output   variable.   Set   equal  to  one  in  C0EF2 

after    integration    of     the     pressure 

th 
distribution  for  the  n    blade  is  completed. 

the  Mach  line  counter.  Incremented  each   time 

a  right-running  Mach  line  is  completed,  it  is 

reset  to  the  value  of  IHAVEP  at   the   leading 

edge  of  a  new  blade. 

an  integer  array.  This  is  the  first  subscript 

of   the   matrix  used  for  storing  the  discrete 

pressure  and  incremental  moment  distributions 

for  each  blade. 

the  value  of  IHAVEP  at  a  surface  or  wake  grid 

point.  When  more  than  one  blade  is  crossed  by 

a  right-running   Mach   line,   LCOUNT   is   set 

equal   to   OLDL   at   the  initial  laft-running 

Mach  line  and   incremented   by   NGHDFN   after 

each  blade. 

the    maximum   value   of    IHAVEP    along    a 

right-running  Mach  line. 

the  value  of  LCOUNT  at  the  last  grid  point  on 

a   blade.  If  IHAVEP  is  greater  than  LIKW,  the 

grid  point  is  in  the  wake. 

identical   to  LB,  this  value  is  passed  to  the 

main  program  via  common. 

an   input   variable.   It   is   the   number   of 

increments  along  a   right-running   Mach   line 

between  blades. 

the  number  of  Ax   increments   along   a   blade 

chord  (of  unit  length) . 

the  total  number  of  grid  points  on   a   blade. 

This   variable   is  passed  to  the  main  program 

via  common. 

an   integer   array   in   COEF  used  to  identify 

successive  grid  points  on  a  blade. 
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NUM 


NUMBLD 
NUMOLD 


NUMPI 
MAXI 


the   number   of  Ax   increments   between   the 

initial  Mach  line  and  the  leading  edge  of  the 

second  blade. 

one  less  than  the  blade  actually  encountered. 

The   program   numbers   blades   from   zero   to 

NUMBLD  -  1. 

the  number  of  blades  in  the  cascade. 

preserves    the    maximum    value    of     NUM 

encountered  along  a  right-running  Mach  line. 

NUMOLD  +  1. 

the   value   of   JLINE  on  a  right-running  Mach 

line  which  intersects  the  leading  edge   of   a 

blade. 


Other  program  variables. 


AI 

A1,A2 

A3,A4 

ANGP12,ANGM12 

ANGP34,ANGM34 

BI 

CR,CI 

DELTA 
DELTAS 

DR,DI 

DSTSTR 
FAZE 

FF1 ,FF2,FF3 


A  ,  the  factor  — kAx 
I  2 

real  and  imaginary  components  of  P 

h 

real  and  imaginary  components  of  P 

phase   angle   between   real   and  imaginary 

components   of   the  lift  force  and  moment  due 

to  plunge 

same  as  for  ANGP12, ANGM 12,  but  for  pitch 

1   M2 

B  ,  the  factor  -k Ax 

I  4  M2-1 

C  ,C  ,  defined  by  Egs.  (IV-13c,d) 
R   I 

the  interblade  phase  angle  in  radians 

the   step   size   along  a  characteristic  (Mach 

line) 

D  ,D  ,  defined  by  Egs.  (IV~13e,f) 
R   I 

the  step  size  along  a  blade 

the  interblade  phase  angle  in  degrees 

the  flutter  frequency  based  on  the   imaginary 

root  and  the  two  real  roots  respectively 
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FS1,FS2,FS3 

FSTRMN 
G1,  ...  G9 

HDSTRL 

K12R,  ...  K56I 

L1#  ...  m 


M1,  ...  M4 


OMEGAA 


the  flutter  speed,  based  on  the  imaginary  root 

and  the  two  real  roots  respectively 

the  input  Mach  number 

variously   defined   in   each  subroutine  which 

solves  for  perturbation  quantities 

one-half  DSTSTR 

defined  in  III-B  for  each  subroutine 

real   and   imaginary   lift   coefficients   for 

plunge  and   pitch   oscillation   respectively. 

Defined  in  Section  IV. 

real  and  imaginary   moment   coefficients  for 

plunge   and   pitch   oscillation  respectively. 

Defined  in  Section  IV 

the  factor  /xr2 


OMEGAH  the  factor  ^-(~) 

PR12,PI12      matrix   storage   for   the   real  and  imaginary 

pressure  distribution  resulting  from  plunge 
PR34,PI34      matrix   storage   for   the   real  and  imaginary 

pressure  distribution  resulting  from  pitch 
P12MAG,P34MAG  absolute   magnitude  of  the  complex  lift  force 

due  to  plunge  and  pitch  rexpectively 
REDFfiQ         the  reduced  frequency 
RINV  (1/k) 

RMR12,RMI12    real   and   imaginary  components  of  the  moment 

distribution  due  to  plunge 
RMR34,RMI34    real   and   imaginary  components  of  the  moment 

distribution  due  to  pitch 
R12MAG, R34MAG   absolute   magnitude   of   the   complex   moment 

resulting  from  plunge  and  pitch  respectively 

RL1,RL2        real  and  imaginary  component  of  L 

h 

RL3rRL4        the  real  and  imaginary  components  of  L 

a 

RM1,RM2        real  and  imaginary  components  of  M 
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RM3,RM4        real  and  imaginary,  components  of  M 

RMUU  the  blade  density  parameter 

RSUBA  the  radius  of  gyration 

S  the  factor  /m*-1 

STGANG         the  input  stagger  angle,  in  degrees 

STEP  amount   by  which  (1/k)  is  incremented  after. a 

call  to  FLUTER 
TNWPST         input  variable  defining  the  vertical  distance 

between  adjacent  blades 
TRNGLH         sin  (a) •DELTAS.    The   distance  between  blades 

divided  by  the  grid  fineness  ratio 
U22R,    . ,  .  C22I  perturbation   quantities  resulting  from  pitch 

at  grid  points  not  on  a  lower  surface 
U33R,  ...  C33I  same   as  above,  but  at  grid  points  on  a  lower 

surface 
U44R,  ...  C44I  perturbation  quantities  resulting  from  plunge 

at  grid  points  not  on  a  lower  surface 
U55R,  ...  C55I  same  as  above,  but  for  grid  points  on  a  lower 

surface 
VRPLNG, VIPLNG   real   and   imaginary  components  of  the  normal 

velocity  perturbation  due  to   plunge   at   the 

leading  edge  of  the  first  blade 
VRPANL, VIPANL  real  and  imaginary  components  of   the   normal 

velocity   perturbation   due   to   pitch  at  the 

leading  edge  of  the  first  blade 

W  the  factor  

M2-1 

WHWA  the  ratio  of  natural  frequencies  of  plunge  to 

pitch 
XI  /ira  {X} 

XN  floating  point  value  of  NUM 

XNEW  blade  position  relative  to  the  leading   edge 

of  the  plade 
XLNGTH         length  of  a  Mach  line  between  adjacent  blades 
XPT  matrix   array   to   store   the   relative  blade 
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position  of  each  grid  point 
XB1,XR2        the  square  roots  of  the  two  real  roots  of  the 

flutter  determinant 
XSUBA  the    position   of   the  center  of   gravity, 

measured  from  the  elastic  axis 
XSUBO  the  elastic  axis  position 
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APPENDIX.  B 


PANEL  FLUTTER 


The   following   variables    appear    in    the   Cylindrical 
Shell/Panel  Flutter  program. 


A2 


BETA 
CPR,CPI 

DELTAS 
DPHIS1 

DPHIS2 

DPHIDX 

DPHIDR 

FACT1 

FACT2 

FSTRMN 
IHAVEP 


the  factor 
the  factor 
the  factor 
the  factor : 


As 

4M 

M 

2/M2-1 

kMAs 


/M2-1 

real  and  imaginary  components  of  the  pressure 

coefficient 

distance  along  a  Mach  line  between  radii 

the  directional  derivative   of  0     along   the 

left-running  Mach  line 

the  directional  derivative   of  0     along   the 

right-running  Mach  line 

0 

x 

r 

n  z 

<¥> 

the  Mach  number 

the   number   of   a  computational  step  along  a 

right-running  Mach  line 

set   equal   to   one   if   the  user  desires  the 

entire  flow   field   printed   out.   Otherwise, 

zero 

consecutive  number  of  the  grid  points  on   the 
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outer   cylinder   surface.   Incremented  by  one 

after  each  call  to  RAD1 

M  the  axial  mode  number,  m 

M1  mode  number  of  the  panel  displacement  used  in 

the  computation  of  Q 

mr 

N  the  circumferential  mode  number,  n 

NGRDFN         the     number     of     increments    along    a 

right-running  Mach  line  between  the  outer  and 

inner  cylinder 

QREAL,QIMAG    real  and  imaginary  components  of  Q 

mr 

REDFRQ         the  reduced  frequency  k 

RI,RO  the  inner  and  outer  cylinder  radii 

TANR  the  flow  tangency  condition 

XPT  axial  position  of  grid  points  on  the  flexible 

panel  surface 
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cccccccccccccccccccccccccccccccceccccccccccccccccccccccccccc 

C  TWO  OEGREE  OF  FREEDOM  PROGRAM 

xccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

THIS  IS  THE  MAIN  PROGRAM 

**#**   NOTE   ***** 

THE  SIZE  OF  THE  ARRAYS  IN  THIS  PROGRAM  ARE  SUBJECT   ' 
TO  THE  FOLLOWING  CONSTRAINT: 

NGRDFN*TAN(A)*(  (  N-l  )  *  (T  AN(  B  i    -  COT(A)  )  +  1/D  ) 
LESS  THAN  TWICE  THE  DIMENSION  OF  THE  ARRAY. 
A   -   MACH  ANGLE 
B   -   STAGGER  ANGLE 
D   -   DISTANCE  BETWEEN  BLADES 
N   -   NUMBER  OF  BLADES 

REAL*8  X,Y,DSTSTR,HDSTRL,TRNGLH,STAG,FSTRMN 

INTtGER  OLDL 

DIMENSION  X(400) ,Y(400) 

DIMENSION  U22R(400),U22I(400),V2  2R(400)  ,  V22K400)  , 
1    C22R(400) ,C22I(400) 

DIMENSION  U33R(4  00) , U33 I ( 400 ) , V33R( 400 ) ,V33I(400J , 
1    C33R ( 400) ,C33 I ( 400) 

DIMENSION  U44R(4  00) ,U44I(400) ,V44R(400) ,  V44K400) , 
1    C44R(400),C44I(400) 

DIMENSION  U5  5R( 400) ,U55 1(400 ),V5  5R (400) ,V55 I (400), 
1    C55R ( 400) ,C55 I ( 400) 

DIMENSION  RL 1(76) ,RL2( 76), RL3( 76), RL4( 76), RMK 76), 

1  RM2( 76) ,RM^( 76) ,RM4( 76) ,ANGP1 2(76) ,ANGM 12(76) , 

2  ANGP34(76),ANGM34(76),P12MAG< 76 ) , P34MAG ( 76 ) , 

3  R12MAG( 76).R34MAG( 76) 
DIMENSION  XPT(4,100) ,PR1 2( 4 , 100) ,PI 12 ( 4, 100) ,PR34( 4, 

COMMON  /BLCK1/TX,Y,REDFRQ,VRPANL,VIPANL,VRPLNG,VIPLNG, 

CCMMGV/BLCK2/  U22R,U2  2I,V2  2R,V2  2I,C2  2R,C22I, 
1    U2RNEW,U2INEWfV2RNEW,V2INEW,C2RNEW,C2INEW 

COMMON  /BLCK3/  U33R ,U33I , V33R, V33I ,C33R , C33 I , 
1    U3RNEW,U3lNEWfV3RNEW,V3INEW,C3RNEW,C3INEW 
1COMMON  /BLCK4/  IHAVEP,NGRDFN,NUM,NUMBLD,NSTPTS,LCOUNT, 

COMMON  /BLCK5/  U44r!u44I , V44R, V44I ,C44R , C44 I , 
I    U4RNEW,U4INEW,V4RNEW,V4INEW,C4RNEW,C4INEW 

COMMON  /BLCK6/  U55R.U55I , V55R , V55 I ,C55R, C55 I , 
1    U5RNEW,U5INEWfV5RNEW,V5INEW,C5RNEW,C5INEW 

C0MM0N/BLCK7/  RL 1 , RL2, RL3 , RL4,RM 1, RM2 , RM3,RM4 , ANGP12 , 
1rnM^SN,1ul,2;J^P^^A^G,^4»P12MAG»Ri2MAG*p34MAG,R34MAG 

COMMON/ BLCK8/  XPT, PR12 , PI 1 2 , PR34 , PI34, I WR I TE, LVEC , 
1    NITOTfNJUNK 

COMMON  /GEOM/  FSTRMN, DSTSTR, HDSTRL ,TRNGLH, STAG, XSUBO, 
1    UtLiA»rAZE 

COMMON/ FLUTR/  RMUU ,XSUBA,RSUBA , WHWA,STEP , I FLUT 

PRINT  OUT  ALL  INPUT  INFORMATION  VIA  INPUT 

WRITE(6,16) 

CALL  INPUT(IHALT) 

f£r-f^!=I  E^AL  1,  TERMINATE.  A  DIMENSION  SIZE  HAS  BEEN 
EXCEEDED. 


IF  (  IHALT  .EQ.  1  )  GO  TO  14 

THE  OUTER  LOOP  IS  USED  FOR  FLUTTER 


ITERATION. 
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1  CONTINUE 
C 

C      INITIALIZE  FLOW  FIELD  QUANTITIES  AND  LOGIC  VARIABLES 
C 

CALL  INITAL 

CALL  COEF 

FAC  =  1 ./DSTSTR 

NEWDST  =  FAC 

LIMIT  =  (  NUMBLD  -  1  l*NSTPTS  +  NEWDST  +  1 

LIMW  =  NEWDST 

MAXI  =  NGRDFN  *  NSTPTS 

NUM  =  0 

NUMOLD  =  0 

NUMPI  =  1 

OLDL  =  0 

ICO  =  0 

IWRITE  =  0 
C 
C      THE  INNER  LOOP  IS  FOR  FLOW  FIELD  CALCULATION 

2  CONTINUE 
C 
C 

c 

C      IF  POINT  IS  ON  THE  INITIAL  RIGHT  RUNNING  MACH  LINE, 

C      CALL  MACHLN. 

C 

IF  (  LCOUNT  .EQ.  0  )  GO  TO  3 
C 

C      IF  THE  POTNT  IS  ON  A  BLADE  OR  IN  A  WAKE,  GO  TO  4 
C 

IF  (  IHAVEP  .EQ.  LCOUNT  )  GO  TO  4 
C 

C      IF  THE  POINT  IS  ON  THE  INITIAL  LEFT  RUNNING  MACH  LINE, 
C      CALL  MACHLN. 
C 

IF  (  IHAVEP  .EQ.  0  )  GO  TO  3 
C 
C      CONDITION  FOR  THE  FIRST  ENCOUNTER  WITH  A  BLADE 


CALL  COMPXY 


C 

IF  ((  IHAVEP  .EQ.  NUMPI-NSTPTS  )  .AND. 
1    (  JLINE  .EQ.  MAXI  ))  GO  TO  10 
C 

CALL  GENFPT 

IF  (  IHAVEP  .EG.  LIMIT  )  GO  TO  7 

IHAVEP  =  IHAVEP  +  1 

GO  TO  2 
C 

3  CALL  MACHLN 

IF  (  IHAVEP  .EQ.  LIMIT  I  GO  TO  7 
IHAVEP  =  IHAVEP  *  1 
GO  TO  2 
C 

4  ICC  =  1 
C 

C      IF  THE  POINT  IS  IN  THE  WAKE,  CALL  WAKE 

IF  (  IHAVEP  .GT.  LIMW  )  GO  TO  5 
CALL  TOP 
CALL  BOTTOM 
CALL  C0EF1 
IHAVEP  =  IHAVEP  +  1 
IF  (  NUM  .EQ.  0  I  GO  TO  2 
NUM  =  NUM  -  1 
LCOUNT  =  LCOUNT  ♦  NGRDFN 
LIMW  =  LIMW  -  NSTPTS 
GO  TO  2 
C 

5  CALL  WAKE 
CALL  C0LF2 

C 

C  IF    IWRITE    EQUAL       1       ON    RETURN    FROM    C0EF2,     IS    THE 
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C      PRESSURE  DISTRIBUTION  FOR  THIS  BLADE  DESIRED? 
C 

IF  (  IWRITE  .NE.  1  I  GO  TO  6 

M  =  NUM  +  1 

IF  (  (  M  .EQ.  IPRESi  )  .OR.  (  M  .EQ.  IPRES2  I  ) 
1  GO  TO  14 

6  IWRITE  =  0 

IF  (  IHAVEP  .EQ.  LIMIT  )  ICO  =  0 
C 

C      IF  THIS  IS  THE  LAST  POINT  TO  BE  COMPUTED  IN  THE  FLOW 
C      FIELD  AREA,  TERMINATE. 
C 

IF  ((  NUM  .GE.  NUMBLD-1  ) . AND. ( I HAVEP  .EQ.  LIMIT  )) 
I  GO  TO  11 

C 

C      IF  THIS  IS  THE  LAST  GRID  POINT  ON  A  RIGHT-RUNNING 
C      MACH  LINE,  SWITCH  THE  THE  LOGIC  VARIABLES  AND  GO  TO 
C      THE  NEXT  MACH  LINE.  OTHERWISE,  CONTINUE. 
C 

IF  (  IHAVEP  .EQ.  LIMIT  )  GO  TO  7 

IHAVEP  =  IHAVEP  +  1 

IF  (  NUM  .EQ.  0  )  GO  TO  2 

LCOUNT  =  LCOUNT  <-  NGRDFN 

NUM  =  NUM  -  1 

GO  TO  2 
C 

7  CONTINUE 
C 

C      IS  THIS  THE  INITIAL  RIGHT-RUNNING  MACH  LINE? 
C 

IF  (  LCOUNT  .EQ.  0  )  GO  TO  8 
C 

C      THIS  IS  THE  LAST  GRID  POINT  ON  A  RIGHT-RUNNING  MACH 
C      LINE.  STORE  THE  PERTURBATION  QUANTITIES  JUST  COMPUTED 
C      AND  GO  TO  THE  NEXT  MACH  LINE. 
C 

I  =  IHAVEP  ♦  1 

U22R(I)  =  U2RNEW 

U22KI)  =  U2INEW 

V22RU1  =  V2RNEW 

V22KI)  =  V2INEW 

C22RU)  =  C2RNEW 

C22KI)  =  C2INEW 
C 

U44R(I)  =  U4RNEW 

U44HIJ     =    U4INEW 

V44R(I)     =    V4RNEW 

V44KI)     =    V4INEW 

C44R( I)     =    C4RNEW 

C44KI)  =  C4INEW 
C 

C      SWITCH  THE  LOGIC  VARIABLES  TO  START  THE  COMPUTATION. 
C 

8  IHAVEP  =  0 
NUM  =  NUMOLD 
ICO  =  0 

JLINE  =  JLINE  +  1 
C 

C      IF  THIS  IS  NOT  THE  FIRST  BLADE,  RESET  THE  LOGIC  SWITCH 
C      VARIABLES  LCOUNT,  LIMW,  OLDL. 


C 


IF  {  NUMOLD  .GT.  0  J  GO  TO  9 
LCCUNT  =  LCOUNT  +    1 
GO  TO  2 

9  LCOUNT  =  OLDL 

LIMW  =  NUMOLD*NSTPTS  +  NEWDST 
OLDL  =  OLDL  y    1 
GO  TO  2 

10  CALL  NEWBLD 
ICO  =  1 
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CALL  C0EF1 

NUMOLD  =  NUM 

LCOUNT  =  IHAVEP  +  NGRDFN 

JLINE  =  IHAVEP 

IHAVEP  =  IHAVEP  ♦  1 

OLDL  =  IHAVEP 

NUMPI  =  NUM  ♦  1 

MAXI  =  NGRDFN  *  NUMPI*NSTPTS 

NUM  =  NUM  -  1 

GG  TO  2 

11  CGNTINUE 

IF  (  IFLUT  .EQ.  1  )  GO  TO  14 

CALL  FLUTER(REDFRQ,NUMBLD,ITO, I  HALT) 

IF  (  IHALT  .EQ.  1  )  GO  TO  14 

RINV  =  1. /REDFRQ 

IF  (  ITO  .EQ.  1  )  GO  TO  13 

RINV    =    RINV    +    STEP 

REOFRQ    =    l./RINV 

GO    TO    1 

12  RINV  =  RINV  -  STEP 
STEP  =  STEP/5. 
WRITE(6,17)  STEP 
RINV  =  RINV  +  STEP 
REDFRQ  =  l./RINV 
GO  TO  1 

13  CONTINUE 

14  CONTINUE 
WRITE(6,18) 

DO  15  N  =  ltNUMBLD 

15  WRITE(6,19)    N,P12MAG(N) ,ANGP12(N),R12MAG(N)  , 

1  ANGM12(N) ,P34MAG(N) ,ANGP34(N) ,R34MAG(N) , 

2  ANGM341N) 
N    =    14 

WRITE (6, 100) FAZE, RL1(N),RL2(N/,RM1(N),RM2(N),RL3(N), 
1  RL4(N> ,RH3(N) ,RM4(N) 

100    FORMAK »Oe ,5X,9F12.4) 
WRITE(6,16) 
STOP 

16  FORMAT (• 1« ) 

17  FGRMAT( »0* ,35X,« ***   NEW  STEP  =  «,F6.4,'   ***«) 

18  FORMAT( »0« ,5X,« BLADE     P12MAG     ANGP12     R12MAG», 

1  •     ANGM12     P34MAG     ANGP34     R34MAG     8, 

2  •ANGM34',//) 
19F0RMAT(«  *f4XfI5f lX,8F10.4f/J 

END 


C  THIS  IS  SUBROUTINE  FLUTER 

C 

SUBROUTINE  FLUTER (REDFRQ , NUMBLD, ITO, IHALT) 

REAL  L1,L2,L3,L4,M1,M2,M3, M4 

DIMENSION  RL 1(76) ,RL2( 76), RL3( 76), RL4( 76), R Ml (76), 

1  RM2( 76) ,RM3(76) ,RM4(76) , ANGP 12( 76 ) , ANGM1 2( 76) , 

2  ANGP 34 ( 76) ,ANGM34(76) ,P12MAG( 76) ,P34MAG( 76) f 

3  R12MAG( 76),k34MAG(76) 

COMMON/ BLCK7/  RL 1 , RL2 , RL3 , RL4, RM 1, RM2, RM3 , RM4 , ANGPi 2 , 
1    ANGM12tANGP34,ANGM34,P12MAG,R12MAG,P3  4MAG,R34MAG 
COMMON/ FLUTR/  KMUU  ,X SUB A, RSU6A , WhWA, STE P , I FLUT 
DATA  IPRINT  /l/ 
C 

IF  (  IPRINT  .EQ.  2  )  GO  TO  1 
IPRINT  =    2 
WRITE(6,115) 
WRITE(6,U6) 
I  CONTINUE 


I  TO  =  0 

II  =  NUMBLD 
C 

C      FACTOR  IS  USED  TO  MAKE  LIFT  AND  MOMENT  VALUES  AGREE 
C      WITH  GARRICK  AND  RUBINOW  NOTATION. 


C 


C 


FACTOR  =  0.5/(REDFRQ*REDFRQ) 


LI  =  RLK  II)*FACTOR 

L2  =  RL2( I I)*FACTOR 

L3  =  RL3( II)*2.*FACTQR 

L4  =  RL4(II )*2.*FACTOR 

Ml  =  RMK  I  I)*2.*FACTOR 

M2  =  RM2( I I)*2.*FACTOR 

M3  =  RM3(II)*<V.*FACTOR 

M4  =  RM4( II J *4.*F ACTOR 
C 
C      DEFINE  CONSTANTS  TO  BE  USED  IN  THE  FLUTER  SOLUTION. 


RINV  =  l./REDFRQ 
W2  =  WHWA*WHWA 
R2  =  RSUBA*RSUBA 
X2  =  XSUBA*XSU8A 
OMEGAA  =  RMUU*R2 
OMEGAH  =  RMUU*W2 


c 

DI    =    L1*M4    - 

-    L4*M1    +    L2*M3    -    L3*M2 

c 

DR   =    L1*M3    - 

-    L3*M1    -    L2*M4    +    M2*L4 

CR  =  RMUU*(  XSUBA*(  Ml  +  L3  )  -  (  M3  -  OMEGAA  )  - 
1    (  L1*R2  +  RMUU*X2  )  J  +  DR 

CI  =  RMUU*<  XSUSA*<  M2  f  L4  I  -  M4  -  L2*R2  )  +  DI 
C 

IF  (  WHWA  .EQ.  0.  )  GO  TO  20 
C 

C      DEFINE  CONSTANTS  FOR  QUADRATIC  EQUATION  OF  REAL 
C      COMPONENT  AND  FOR  THE  IMAGINARY  COMPONENT. 
C 

B  =  (LI  -  RMUU  ) /OMEGAH  +  M3/0MEGAA  -  1. 

C  =  CR/(0MEGAA*0MEGAH1 

D  =  L2* OMEGAA  +  0MEGAH*M4 
C 

C      SOLVE  FOR  AND  TEST  THE  DISCRIMINANT 
C 

DISCRM  =  B*B  -  4.*C 

IF  (  DISCRM  .LT.  0.  )  GO  TO  40 
C 

C      SOLUTION  OF  IMAGINARY  COMPONENT 
C 

XII  =  -  CI/D 

IF  (  XII  .LT.  0.  )  GO  TO  41 

XI  =  SQRT(XII) 
C 

C      SOLUTION  OF  REAL  COMPONENT 
C 

XR11  =  C  -  B  +  SQRT(DISCRM)  )*0.5 

XR21  =  (  -  8  -  SQRT(DISCRM)  )*0.5 

IF  (  (  XR11  .LT.  0.  )  .OR.  (  XR21  .LT.  0.  )  )  GO  TO  42 

XR1  =  SQRT(XRll) 

XR2  =  SQRTUR21) 
C 

C      CALCULATION  OF  FLUTTER  FREQUENCY  ( FF )  AND  SPEED  (FSi 
C      BASED  ON  ALL  THREE  ROOTS. 
C 

FF1  =  l./XI 

FS1  =  FF1*RINV 

FF2  =  1./XR1 

FS2  -  FF2*RINV 

FF3  =  1 ./XR2 

FS3  -  FF3*RINV 

WRITE (6, 100 i     REDFRQtRINV,XI ,XR 1  ,XR2 ,FF 1 , FS1 ,FF2,FS2, 
1    FF3,FS3 
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c 

C      SOLUTION  WHEN  RATIO  OF  NATURAL  FREQUENCIES  IS  ZERO 


20  F  =  0MEGAA*(L1  -  RMUU) 
C 
C      SOLUTION  FOR  IMAGINARY  ROOT 


XII  =  -  CI/(OMEGAA*L2) 

IF  (  XII  .LT.  0.  J  GO  TO  41 

XI  =  SQRTIXII) 
C 

C      SOLUTION  FOR  REAL  ROOT 
C 

XR11  =  -  CR/F 

IF  (  XR11  ,LT.  0.  J  GO  TO  42 

XR1  =  SQRT(XRll) 
C 

C      CALCULATION  OF  FLUTTER  FREQUENCY  ( FF )  AND  SPEED  (FS) 
C      BASED  ON  THE  REAL  AND  IMAGINARY  COMPONENTS. 


C 


FF1  =  l./XI 

FSl  =  FF1*RINV 

FF2  =  1./XR1 

FS2  =  FF2*RINV 

WRITE (6, 117)  REDFRQ,RINV,XI ,XR 1 , FF1 ,FS 1 , FF2 ,FS2 

RETURN 

40  WRITE(6,110)  REDFRQ,RINV 
I  TO  =  1 

RETURN 

41  WRITE(6,111J  REDFRQ,RINV 
IHALT  =  1 

RETURN 

42  WRITE(6,112)  REDFRQ,RINV 
IHALT  =  1 

RETURN 
100  FORMAT* '0' ,5X,5F10. 4, •    **«,6F10.4) 

110  FORMATMO',  5X  ,  2F1  0.4,  5X  , •  DI  SC  RM  <  0«) 

111  FORhATCO',  5X.2F10.4, 15X,'XIMAG  <  0",//) 

112  FORMATCO",  5X , 2F 1 0. 4, 15X, • XRE AL  <  0',//) 

115  FORMAT! '  0« ,20X,« FLUTTER  POINT  IS  THE  INTERSECTION1, 

1  •  OF  THE  REAL  AND  IMAGINARY  »t/»  •, 20X, "CURVE S  •, 

2  •  OF  SQRT(X)  WHEN  PLOTTED  VERSUS  il/K).   AT  THIS  ', 

3  •POINT:',/'  '  ,30X,  'FLUTTER  FREQUENCY  ( FF )  =  «? 

4  '1/SQRTiX) * ,/•  ', 'FLUTTER  SPEED  (FS)  =  FF/K',/) 

116  FORMATC  • , 30X, • SQRT «  , 5X, • SQRT' , 6X, '  SQRT • , 1 2X, ' BAS ED' , 

1  •  ON  XIMAG'  ,5Xf  '  BASED  ON  XRE  ALU) ',  3X  ,' BASED  ON  e, 

2  'XREAL(-) ',/'  ' ,11X, »K' ,8X,«1/K»,7X, • (XI )« ,5X, • (XR+ 

3  )  •  ,5X,» (XR-) • ,12X,°FF' ,8X,' FS' ,8X, 'FF« ,8X,'FS8 ,8X, 

4  'FF' ,8X,'FS' ,//) 

117  FORMAT! «0« ,5X,4F10. 4, 10X,«    **',4F10.4) 
END 


C  THIS  IS  SUBROUTINE  INPUT 

C 

SUBROUTINE  INPUT(IHALT) 
C 

C      SUBROUTINE  INPUT  READS  ALL  GEOMETRIC  AND  AERODYNAMIC 
C      INPUT  DATA. 
C 


REAL*8  X,Y,DSTSTR,HDSTRL,TRNGLH, ST  AG, FST RMN , FMANGL , 

1  TNW DST, DE LTA S, X LNGTH,STGR, SO LI D,STGANG , ARG, FACTOR, 

2  ANG 

DIMENSION  X(400) ,Y(400) 

COMMON  /BLCK1/  X,Y ,REDFRQ, VRPANL » VIPANL . , RP LNG, VI PLNG , 
1    S  ,  to  ,  a  I  ,  B  I 

COMMON  /BLCK4/  I HAV EP, NGRDFN, NUM, NUMBLD , NSTPTS , LCOUNT, 
1    JLINE,IC0,IPRES1,IPRES2 
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COMMON  /GEOM/  FSTRMN,DSTSTR ,HDSTRL ,TRNGLHf STAG.XSUBO , 
1    DELTA, FAZE 

COMMON/FLUTR/  RMUU  ,X SUBA , RSUBA , WHWA,STE P, IFLUT 

NAMELIST  /IMAM1/  NGRUF  N  tFSTRMN,  REDFRQ,XSUBO  ,  TNWDST  , 
1    STGANG, FAZE, NUMBLD, I  PRE S 1 , 1 PRES2 

NAMELIST/FLTR/  RMUU , XSU6A , RSUBA , WHWA , STE P , IFLUT 

DATA  MAXDIM/400/,  MAXPTS/1QGV,  MAXBLD/75/,  ICROSS/4/ 
C 

READ(5,NAM1> 

WRITE(6,NAM1 J 

READ(5,FLTR) 

WRITE(6,FLTR) 
C 

FMANGL    =    DARSIN( 1.DG7FSTRMN) 

XLNGTH  =  TNWDST*FSTRMN 

FNGDPT  =  DFLGATI NGRDFN) 

DELTAS  =  XLNGTH/FNGDPT 

TRNGLH  =  TNWDST/FNGDPT 

HDSTRL  =  DELTAS*DCOS(FMANGL) 

DSTSTR  =  HDSTRL  +  HDSTRL 
C 

C      COMPUTATION  OF  THE  COMPATIBLE  STAGGER  ANGLE.  THIS  IS 
C      REQUIRED   BECAUSE  THE  LEADING  EDGE  OF  EACH  BLADE  MUST 
C      BE  ON  A  GRID  POINT. 
C 

ARG  =  STGANG*i.74532925D-02 

STAG  =  TNWDST*DTAN<ARG) 

FACTOR  =  TNWDST/DTAN( FMANGL) 

STGR  =  STAG  -  FACTOR 

NSTPTS  =  STGR/DSTSTR 

IF  (  NSTPTS  .LT.  1  )  NSTPTS  =  1 

STGR  =  DFLOAT(NSTPTS)*DSTSTR 

STAG  =  STGR  +  FACTOR 

ANG  =  DATAN(  STAG/TNWDST  J 

STGANG  =  ANG*57.29577951D0 
C 

SOLID  =  DCOS(ANG)/TNWDST 

S  =  DSQRT{  FSTRMN*FSTRMN  -  l.DO  ) 

W  =  FSTRMN*FSTRMN/(S*S) 

DELTA  =  FAZE*0.i745329E-0l 
C 

WRITE ( 6, 10)NGRDFN,FSTRMN,REDFRQ, TNWDST, SOLI D.XSUBO, 
1    STGANG, DSTSTR, FAZEtNUMBLD 
C 

C      CHECK  WHETHER  OR  NOT  VECTOR  DIMENSIONS  ARE  EXCEEDED. 
C      I8LDPT  IS  THE  TOTAL  NUMBER  OF  GRID  POINTS  ON  A  BLADE, 
C      INTXN  .LE.  (IBLDPT+1)  INDICATES  THAT  A  MACHLINE 
C      CROSSES  TOO  MANY  BLADES,  AND  IPTS  IS  THE  TOTAL  NUMBER 
C      OF  GRID  POINTS  ALONG  THE  MACHLINES.   ICROSS  IS  THE 
C      MAXIMUM  NUMBER  OF  BLADES  THAT  CAN  BE  CROSSED  BY  A 
C      RIGHT-RUNNING  MACH  LINE. 


C 


I  HALT  =  0 

INTXN    =    ICROSS*(     NSTPTS    «•    NGRDFN    ) 

FAC    =    1. /DSTSTR 

IBLDPT    =    FAC    «-    2 

IPTS    =    (    NUMBLD    -    I    )*NSTPTS    +    IBLDPT 

IF    (  IPTS    .GT.    MAXDIM    I    GO    TO    4 

1  IF    (  IBLDPT    .GT.    MAXPTS    )     GO    TO    5 

2  IF    (  NUMBLD    .GT.    MAXBLD    )     GO    TO    6 

3  IF    (  INTXN    .LE.     IBLDPT    +1     J    GO    TO    7 
IF    (  IHALT    .EQ.     1     I     RETURN 

GO    TO    8 

4  WRITE(6f120)  MAXDIM 
IHALT  =  1 

GO  TO  1 

5  WRITE(6,12l)    MAXPTS 
IHALT    =    1 
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GO  TO  2 

WRITE(6,i22)  MAXBLD 
IHALT  =  1 
GO  TO  3 

WRITE(6,123)  ICROSS 
IHALT  =  1 
RETURN 


8  CONTINUE 
RETURN 
10  FORMAT!////1 


1  /'0' ,5X, 'FREEST 

2  'REDUCED  FREQUE 

3  'BETWEEN  AIRFOI 

4  F10.7,/'0' ,5X,« 

5  « AXI S    =     • ,F10.7 

6  »    =     ',F10.7,/'0 

7  'UPPER    AIRLOIL 

8  •  BLADES',/// 
11    FORMATC  1»  ,34X,»  TE 

1  'AMPLITUDES'/// 

2  'UI ' ,T58,« VR"  ,T 

120  FORMAT ( 'O1 ,2X,74( ' 

1  'NUMBER  OF  GRID 

2  ISt/'O' ,2X,74(« 

121  FCRMATCO*  ,2X,74(  • 

1  'NUMBER    OF    POIN 

2  2X,74(«*,)J 

122  FORMAT  (  '0'.,2X,54(» 
1  'NUMBER   OF    BLAD 

123  FORMAT  ( '0' ,2X,77( • 

1  'RIGHT-RUNNING 

2  »   BLADES', /'O' 
END 


',5X,'GRID  FINENESS  INPUT  N 
""'REAM  MACH  NUMBER  =  ',F 
NCY  =  • ,F10. 7t/'0',5X 
LS  =  ' ,F10.7,/"0",5X, 
HORIXONTAL  POSITION  OF 

/•0« ,5Xt 'COMPATABLE  S 
• ,5X, 'DELTA  X  =  ',F10. 
LEADS  LOWER  BY  •  ,F10.4 
) 

IPEL  COMPLEX  PERTURBAT 
•  •  ,1X, 'BLD« ,T10,'X8  ,T 
73, • VI' ,T90, «CR' ,T107, 
*•  ),/"0"  ,5X, "  ERROR  TER 

POINTS  ALONG  A  MACH  L 
*•)  ) 

*■ li/iOt ,5X, 'ERROR  TER 
TS  ALONG  A  BLADE  EXCEE 

)  ,  /•0«  ,5X, 'ERROR  TE 
ES  EXCEEDS  ' ,I4,/«0« ,2 
*"  ),/«0'  ,5X, ' ERROR  TER 
MACH  LINE  CROSSES  MORE 
,2X,77( •*« )) 


UMBER  =  • ,110, 
10. 7, /«0« , 5X, 
•DISTANCE  ', 
SOLIDITY  =  ', 

ELASTIC  ', 
TAGGER  ANGLE', 
7,/'0« ,5X, 
,/'0',I3, 

ION  •  , 

26, 'UR« ,T42, 
•CI" ,/) 
MINATION:  ', 
INE  EXCEEDS  ', 

MINATION:  ', 
DS  ,,I5,/,0"f 

RMINATION:  ', 
X,54( •*« ) I 

MINATION:  ', 
THAN  ' , 12, 


THIS  IS  SUBROUTINE  INITAL 

SUBROUTINE  INITAL 

SUBROUTINE  INITAL  INITIALIZES  ALL  FLOW  FIELD 
QUANTITIES. 


REALMS  X*Y,DSTSTR, 

DIMENSION  X(400) ?Y 

DIMENSION  U22R1400 
1    C22R(400) ,C22I< 

DIMENSION  U33RC400 
1    C33R(400),C33I( 

DIMENSION  U44R(400 
1    C44R(400),C44I( 

DIMENSION  U55R(400 
1    C55RC400) ,C55I< 

COMMON  /BLCKL/  X,Y 
1    S  ,  W ,  A  I  ,  B I 

COMMON  /BLCK2/  U22 
1    U2RNEW,U2INEW, V 

COMMON  /BLCK3/  U33 
1    U3RNEW,U3INEW,V 

COMMON  /BLCK4/  IHA 
1    JLINE,ICO,IPRES 

COMMON  /6LCK5/  U44 
1    U4RNEW,U4INEW,V 

COMMON  /BLCK6/  U55 
1    U5RNEW,U5INEW,V 

COMMON  /GEGM/  FSTRMN 
1    DELTA, FAZE 


AI  =  .5*REDFRQ*SNGLtDSTSTR) 


HDSTRL,TRNGLH, ST AG, FSTRMN 

(400) 

)  ,U22  1(400  ),  V22R(  400)  ,  V22K400)  , 

400) 

),U33I(400),V33R(400) ,V33I(400), 

400) 

),U44I(400),V44R(400) ,V44I(400) , 

400) 

)TU55I(400),V5  5R(400) ,V55I(400) , 

400) 

,REDFRQ,VRPANL,VIPANL,VRPLNG,VIPLNG, 

R,U22I,V2  2R»V2  2I,C22R,C22I, 

2R NE  W , V2 I NE  W , C 2RNE  W , C 2  I NE  W 

R,U3  3I,V33R,V33I,C33R,C33I, 

3KNEW,V3INEW,C3RNEW,C3INEW 

VEP,NGRDFN,NUM,NUMBLD,NSTPTS,LCOUNT, 

1,IPRES2 

R.U44I, V44R, V44 I ,C44 R , C44 I , 

4RNEW,V4INEW,C4RNEW,C4INEW 

R,U5  5-I,V5  5R,  V55I,C55R,C55If 


,DSTSTR,HOSTRL,TR 


LH,STAG,XSUBO, 
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BI    =    .5*AI*W 
X(l)    =    O.ODO 
XI    =    X( 11 
Y(l)    =    O.ODO 


SET    UP    INITIAL    VALUES    FOR    U,V,     AND    C    AT     (0,0) 
PITCH: 


VRPANL  =  -l./S 

VIPANL  =  REDFRQ*XSUBO/S 


PLUNGE: 


VRPLNG  =  0.0 
VIPLNG  =  -REDFRQ/S 


INITIAL  VALUES  OF  PERTURBATION  QUANTITIES   RESULTING 
FROM  PITCH  OSCILLATION  AT  THE  LEADING  EDGE  OF  THE 
FIRST  BLADE;  UPPER  SURFACE  (22),  LOWER  SURFACE  (33) 


U22R(1 
U22K1 
V22RU 
V22K  1 
C22R(1 
C22K1 
U33R( L 
U33K1 
V33RU 
V33K1 
C33R( 1 
C33K  1 


-VRPANL 

-VIPANL 

-U22RU) 

-U22I ( 1) 

-U22R(  I) 

-U22K  1) 

VRPANL 

VIPANL 

VRPANL 

VIPANL 

-VRPANL 

-VIPANL 


INITIAL  VALUES  OF  PERTURBATION  QUANTITIES   RESULTING 
FROM  PLUNGE  OSCILLATION  AT  THE  LEADING  EDGE  OF  THE 
FIRST  BLADE;  UPPER  SURFACE  (44),  LOWER  SURFACE  (55) 


U44RC1) 
U44K1) 

V44R( 1) 
V44K  1) 
C44R( 1) 
C44K  1) 
U55R(1) 
U55K1) 
V55R(1) 
V55K1  ) 
C55R( 1) 
C55K  1) 


-VRPLNG 
-VIPLNG 
-U44R( I) 
-U44K  1) 
-U44R( 1) 
-U44K  1) 
VRPLNG 
VIPLNG 
VRPLNG 
VIPLNG 
-VRPLNG 
-VIPLNG 


JLINE  =  0 
IHAVEP  *    1 
LCOUNT  =  0 
RETURN 
100  F0RMATI2H  I,7E16.7,/) 
END 


THIS  IS  SUBROUTINE  COMPXY 
SUBROUTINE  COMPXY 
COMPXY  COMPUTES  THE  (X,Y)  GRID  POINT  POSITION 

REAL*8  X,Y,DSTSTR,HDSTRL,TRNGLH,STAG,FSTRMN 

DIMENSION  X(400) , Y(400) 

CCMMON  /BLCK1/  X , Y , REDFRQ, VRPANL ,V IP ANL , VRPLNG, VI  PLNG, 
1  S ,  W ,  A  I ,  B  I 

COMMON  /BLCK4/  IHAVEP, NGROFN, NUM, NUMBLD, uST PTS , LCOUNT, 
1  JLINE, ICO, IPRESi,IPRES2 
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CCMMON  /GEOM/  FS  FRMN,  DSTST  R,  HOS  TRL  ,TRNGLH,  STAG,  XSUBO, 
1    DELTA, FAZE 

I  =  IHAVEP  *•  1 

IF  GRID  POINT  IS  ON  THE  INITIAL  LEFT  RUNNING  MACH 
LINE,  GO  TO  1. 

IF  ((  IHAVEP  .EQ.  0  J.AND.l  LCOUNT  .NE.  0  ) )  GO  TO  1 
X(I)  =  X(  IHAVEP)  *■    HDSTRL 
Yd)  =  Y(  IHAVEP)  -  TRNGLH 
RETURN 

1  X( I)  =  X(l)  +  HDSTRL 
Yd)  =  Y(l)  ♦  TRNGLH 
RETURN 
END 


THIS  IS  SUBROUTINE  MACHLN 

SUBROUTINE  MACHLN 

MACHLN  COMPUTES  THE  VALUE  OF  U,V,  AND  C  ALONG  THE 
INITIAL  MACH  LINE  AT  THE  GIVEN  GRID  POINT. 


REAL*8 

DIMENSI 

D I  ME  NS I 
1    C22R 

DIMENSI 
1    C33R 

DIMENSI 
1    C44R 

DIMENSI 
1    C55R 

CCMMON 
1   S ,  w 

COMMON 
1    U2RN 

COMMON 
1    U3RN 

COMMON 
1    JLIN 

COMMON 
1    U4RN 

COMMON 
1    U5RN 


X,Y 

ON  X(40 
ON  U22R 
(400) ,C 
ON  U33R 
(400) ,C 
ON  U44R 
(400) ,C 
ON  U55R 
(400) ,C 
/BLCK1/ 
AI  ,31 
/BLCK2/ 
EW,U2IN 
/BLCK3/ 
EW,U3IN 
/BLCK4/ 
E,  ICO, I 
/3LCK5/ 
EW,U4IN 
/BLCK6/ 
EW,U5IN 

1 


0) ,Y(400) 

(400)  ,U 22  I (400),  V22R{400)  ,  V22K400)  , 
22  1(400) 

(400) ,U33I(400),V33R(400) ,V33I(400) , 
.33  1(400) 

(400)  ,U44I(400),  V44RI400)  ,  V44K400)  , 
.44  1(400) 

(4  00),U55K400)«V55R(400)  ,  V55K400)  , 
551(400) 
X,Y,REDFRQ,VRPANL,VIPANL,VRPLNG,VIPLNG, 

U2  2R,U22I,V2  2R,V22I,C22R,C22I, 
EW , V2RNEW , V2 I  NEW , C2RNEW, C2 INEW 

U33R,U33I,V33Ri V33I ,C33R , C33I » 
EWt V3RNEW,V3 INEW fC3RNEW,C3 INEW 

IHAVEP,NGRDFN,NUM,NUMBLD»NSTPTS,LCOUNT, 
PRES1,IPRES2 

U44R,U44[ , V44R, V44 I ,C44R , C44I , 
EW , V4RNE  W , V4 1  NEW , C4RNE  W , C  4 1 NE W 

U5  5R,U55I, V5  5fc,V5  5I,C55kf C55If 
EW,V5R.MEW,V5INEW,C5RNEW,C5INEW 


I  =  IHAVEP 
XI  =  X(I) 

ARG  =  REDFRQ*W*XI 
U  =  COS(ARG) 
V  =  SIMARG) 

IF  POINT  IS  ON  THE  INITIAL  RIGHT  RUNNING  MACH  LINE, 
GO  TO  2 

IF  (  IHAVEP  .NE.  0  )  GO  TO  2 

THIS  PORTION  IS  FOR  PITCH 


VIPANL#V 
VRPANL*V 


U2RNEW  =  -VRPANL*U 
U2INEW  =  -VI PANL*U 
V2RNEW  =  -U2RNEW 
V2INEW  =  -U2INEW 
C2RNEW  =  -U2RNEW 
C2INEW  =  -U2INEW 


THIS  PORTION  IS  FOR  PLUNGE 
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-VRPLNG*U 
-VIPLNG*U 

-U4RNEW 
-U4INEW 
-U4RNEW 
-U4INEW 


VIPLNG*V 
VRPLNG*V 


U4RNEW 
U4INEW 
V4RNEW 
V4INEW 
C4RNEW 
C4INEW 
RETURN 


THIS    PORTION    IS    FOR    PITCH 


U22RU) 
U22KI) 
V22RU ) 
V22KI) 
C22R( I) 
C22H  I) 


VRPANL*U 
VIPANL*U 
U22RU  ) 
U22K  I  ) 
-U22RC  I) 
-U22KI  ) 


VIPANL*V 
VRPANL*V 


THIS    PORTION    IS    FOR    PLUNGE 


U44RU1  = 

U44KI)  = 

V44R(I)  = 

V44KI)  = 

C44RU)  = 

C44KI)  = 
RETURN 
100    FORMAT (2H 
END 


VRPLNG*U 
VIPLNG*U 
U44R( I ) 
U44K  I  ) 
-U44R( I ) 
-U44K  I) 


M,7E16.7,/) 


+  VIPLNG*V 
-  VRPLNG*V 


THIS  IS  SUBROUTINE  TOP 
SUBROUTINE  TOP 
TOP  COMPUTES  U,V  AND  C  AT  AN  UPPER  SURFACE  POINT 

REAL*8  X,Y,DSTSTR,HDSTRL,TRNGLH,STAG,FSTRMN 

REAL  K12R,K12I ,K34R, K34I f K56R, K56I 

DIMENSION  X<400J  ,Y(400) 

DIMENSION  U22RC400) , U22I ( 400 ) , V2  2R(400) ,  V22H400) , 
1    C22R(400) ,0221(400) 

DIMENSION  U44R(400)  ,U44I  (  400  )  ,  Vf4R(400  >  ,  V44K400)  , 
1    C44R(400) ,C44I<400) 

COMMON  /BLCKI/  X  ,  Y  ,REDFRQ, VRPANL , VI PANL , VRP LNGj VI PLNG, 
1    S,W,AI,BI 

COMMON  /BLCK2/  U22R,U22I , V2 2R, V22I ,C22R, C22 I , 
1    U2RNEW,U2INEW,V2RNEW,V2INEW,C2RNEW,C2INEW 

COMMON  /BLCK4/  I HAVEP, NGRDFN.NUM, NUMBLD , NSTPTS ,LCOUNT, 
1    JUNE,  ICO,  IPRESli  IPRES2 

COMMON  /BLCK5/  U44R, U44I , V44R, V44I ,C44R , C44I , 
1    U4RNEW,U4INEW,V4RNEW,V4INEW,C4RNEW,C4INEW 

COMMON  /GEOM/  FS TRMN , DSTSTR , HDSTRL ,TRNGLH,STAG, XSU80 , 
1    DELTA, FAZE 

XN  =  FLOAT(NUM) 

J  =  IH4VEP 

I  =  I HAVEP  +  1 

XNEW  =  SNGL(XU))  -  SNGL( STAG) *XN 

ARG  =  DELTA*XN 

THIS    PORTION    IS    FOR    PITCH 

DEFINE    CONSTANTS.     (J)     SUBSCRIPT    REFERS    TO    PERTURBATION 
QUANTITIES    ON    THE    PREVIOUS    MACH    LINE    AT    GRID    POINT 
IHAVEP.     (i\IEW)     INDICATES    QUANTITIES    ON    THE    CURRENT    MACH 
LINE    AT     IHAVEP 

KI2R  =  U22RU)  +  C22R(J)  <-  AI*U22MJ) 
KI2I  =  U22HJ)  +  C22KJ)  -  AI*U22R(JJ 
K34R   =    -(    COS(ARG)    -    (     XNEW    -    XSUBO    )*REDFRQ* 
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1  SIN(ARG)  )/S 

K34I  =  -(  SIN(ARG)  *    (  XNEW  -  XSUBO  )*REDFRQ* 

1  COS(ARG)  )/S 

K56R  =  U2RNEW  ♦  V2RNEW  ♦  BI*<  U2INEW  -  C2INEW  ) 

K56I  =  U2INEW  +  V2INEW  -  BI*(  U2RNEW  -  C2RNEW  ) 

Gl  =  1.  -  AI*BI 

G2  =  BI  ♦  BI 

G3  =  G1*G1  f  G2*G2 

G4  =  K56R  -  K34R  -  BI*K12I 

G5  =  K56I  -  K34I  +  BI*K12R 

NOW  STORE  THE  PERTURBATION  QUANTITIES   FOR  GRID  POINT 
IHAVEP  ON  THIS  MACH  LINE 


U22R(J) 
U22K  J) 


=  U2RNEW 
=  U2INEW 


V22RUI  =  V2RNEW 

V22I (JJ  =  V2INEW 

C22R(J)  =  C2RNEW 

C22KJ)  =  C2INEW 

COMPUTE  PERTURBATION  QUANTITIES  AT  THE  CURRENT  GRID 
POINT  (IHAVEP+1),  AND  TEMPORARILY  STORE  AS  (NE^J 

U2RNEW  =  (  G4*G1  *-  G5*G2  )/G3 

U2INEW  =  (-G4*G2  *-  G5*G1  >/G3 

V2RNEW  =  K34R 

V2INEW  =  K34I 

C2RNEW  =  K12R  -  U2RNEW  ♦  AI*U2INEW 

C2INEW  =  K12I  -  U2INEW  -  AI*U2RNEW 

THIS  PORTION  IS  FOR  PLUNGE.  STEPS  ARE  THE  SAME  AS  FOR 
PITCH 

K12R    =    U44RU)     +    C44R(J)     *•    AI*U44l(J) 

KL2I    =    U44KJ)     +    C44KJ)     -    AI*U44R(J) 

K34R    =    REDFRQ*SIN( ARG)/S 

K34I    =    -REDFRQ*COS(ARG}/S 

K56R    =    U4RNEW    +    V4RNEW    +    BI*(     U4INEW    -    C4INEW    ) 

K56I    =    U4INEW    +    V4INEW    -    BI*{     U4RNEW    -    C4RNEW    ) 

G4   =    K56R    -    K34R    -    BI*K12I 

G5  =  K56I  -  K34I  +  BI*K12R 

U44RU)  =  U4RNEW 

U44KJ)  =  U4INEW 

V44RU)  =  V4RNEW 

V44KJ)  =  V4INEW 

C44R1J)  =  C4RNEW 

C44HJ)  =  C4INEW 


U4RNEW 
U4INEW 

V4RNEW 
V4INEW 
C4RNEW 

C4INEW 


(  G4*G1 

(-G4*G2 

K34R 

K34I 

K12R  -  U4RNEW 

K12I  -  U4INEW 


G5*G2  J/G3 
G5*G1  )/G3 


AI*U4INEW 
AI*U4RNEW 


N  =  NUM  +  1 
RETURN 
100  F0RhAT(2H  T , I  2 , 7E 16.7, / ) 
END 


THIS  IS  SUBROUTINE  BOTTOM 
SUBROUTINE  BOTTOM 
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REAL*8 

REAL    Kl 

DIMENSI 

D I  MENS  I 
1  C22R 

DIMENSI 
1         C33R 

0 1  ME  NS  I 
1  C44R 

DIMENSI 
1         C55R 

COMMON 
1  S,W 

COMMON 
1  U2RN 

COMMON 
1  U3RN 

COMMON 
1  JLIN 

COMMON 
1  U4RN 

COMMON 
i         U5RN 


X,Y 

2R,K12I 
GN    X(400 
ON    U22R( 
(400) ,C2 
ON    U33R( 
(400) ,C3 
ON    U44R( 
(400) ,C4 
ON    U55R( 
(400) ,C5 
/BLCKi/ 
AI  ,BI 
/BLCK2/ 
EWtU2INE 
/BLCK3/ 
EW,U3INE 
/BLCK4/ 
E,ICO,IP 
/BLCK5/ 
EW,U4INE 
/BLCK6/ 
EW,U5INE 


K34R,K34I,K56R,K56I 
)  ,Y(400) 

400)  ,U22I(400),V22R(400)  ,  V22K400)  , 
2  1(400) 

4  00)  ,  U33K400),  V33R(400)  ,  V33K400)  , 
31(400) 

400),U44I(400),  V44R(400)  ,  V44K400)  , 
41(400) 

400)  ,U55I(400)fV55R(400)  .V55K400I  t 
51(400) 
,Y, REDFRG,VRPANL,VIPANL,VPPLNG,VIPLNG, 

U22R,U22I,V2  2R,V22I,C22R,C22I, 

W, V2RNEW,V2INEW,C2RNEW,C2INEW 

U3  3R,U33I,V3  3R, V33  I  ,C33R , C33 I , 

W,V3RNEW,V3INEW,C3RNEW,C3INEW 

IHAVEP,NGRDFN,NUM,NUM8LD,NSTPTS,LC0UNT, 

RESlt IPRES2 

U44R,U44I, V44R, V44 I ,C44R , C44 I , 

W,V4RNEW,V4lNEw* C4RNEW,C4INEW 

U55R,U5  5I,V5  5R, V55I,C55RtC55i, 

W,V5RNEW,V5INEW,C5RNEW,C5INEW 


BOTTOM    COMPUTES    U,V    AND    C    AT    A    LOWER    SURFACE    POINT 

J    =    IHAVEP 

I    =    IHAVEP    4-    1 

xi  =  xm 

THIS    PORTION    IS    FOR    PITCH 

DEFINE    CONSTANTS.     (J)     SUBSCRIPT    REFERS    TO    PERTURBATION 
QUANTITIES    ON    THE    PREVIOUS    MACH    LINE    AT    GRID    POINT 
IHAVEP.     (NEW)     INDICATES    QUANTITIES    ON    THE    CURRENT    MACH 
LINE    AT    IHAVEP 


K12R  =  U33R(J) 
K12I  =  U33KJ) 
K34R  =  V2RNEW 
K34I  ~  V2INEW 
K56R  =  U22R( I) 
K56I    =    U22K  I) 


C33R( J) 
C33K  J) 


V22R( I) 
V22I  (  I) 


i-    AI*U33I(J) 
-AI*U33R( J) 


BI*(     U22KI)     -    C22KIJ     ) 
BI*(     U22R(I)     -   C22RU)     ) 


Gl  =  1.  -  AI*BI 

G2  =  BI  ♦  BI 

G3  =  G1*G1  ♦  G2*G2 

G4  =  K56R  +    K34R  -  BI*K12I 

G5  =  K56I  +  K34I  +    BI*K12R 

COMPUTE  PERTURBATION  QUANTITIES  AT  THE  CURRENT  GRID 
POINT  (IHAVEP+1),  AND  TEMPORARILY  STORE  AS  (NEW) 

U3RNEW  =  (  G4*G1  ♦  G5*G2  ) /G3 

U3INEW  =  (~G4*G2  +  G5*G1  )/G3 

V3RNEW  =  K34R 

V3INEW  =  K34I 

C3RNEW  =  K12R  -  U3RNEW  *•  AI*U3INEW 

C3INEW  =  Ki2I  -  U3INEW  -  AI*U3RNEW 

STEPS  ARE  THE  SAME  AS  FOR 


AI*Uf55I(  J) 
-AI*U55R( J) 


c 
c 
c 
c 

c 

THIS  PORTION  IS  FOR  PLUNGE 
PITCH 

K12R  =  U55R(J)  *■    C55R(J)  *■ 

K12I  =  U55KJ)  +  C55HJ)  - 

K34R  =  V4RNEW 

K34I  =  V4INEW 

K56k  ^  U44K(  I  )  -  V44R<  I  )  + 

K56I  =  U44H  I  )  -  V44I  {  I  )  - 

G4  =  K56R  +  K34R  -  BI*K12I 

8I*<     U44IU) 
BI*{     U44RCI) 


C44K  I  ) 
C44R(I  ) 
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G5  =  K56I  ♦  K34I  +  BI*K12R 


U5RNEW 

U51NEW 
V5RNEW 
V5INEW 
C5RNEW 
C5INEW 
RETURN 

100  FGRMAT(2H  B,7E16.7,/) 
END 


(    G4*G1 

(-G4*G2 

K34R 

K34I 

K12R    -    U5RNEW 

K12I    -    U5INEW 


G5*G2    )/G3 
G5*G1    J/G3 


AI*U5INEW 
AI*U5RNEW 


THIS    IS    SUBROUTINE    GENFPT 
SUBROUTINE    GENFPT 
GENFPT    COMPUTES    U,V,     AND    C    AT    A    GENERAL     FIELD    POINT 


REAL*8 
REAL    Kl 
DIMENSI 
D I  MENS  I 

1    C22R 
DIMENSI 

1    C33R 
DIMENSI 

1    C44R 
DIMENSI 

1    C55R 
COMMON 

1  s.w 

COMMON 
1    U2RN 

COMMON 
I    U3RN 

COMMON 
1    JLIN 

COMMON 
1    U4RN 

COMMON 
1    U5RN 


X,Y 

2R,K12I 
ON  X(40 
ON  U22R 
(400) ,C 
ON  U33R 
<400) ,C 
ON  U44R 
(400) ,C 
ON  U5  5R 
(400) ,C 
/DLCK1/ 
AI  ,BI 
/6LCK2/ 
EW,U2IN 
/BLCK3/ 
EWf U31N 
/BLCK4/ 
E,  ICO,  I 
/BLCK5/ 
EWfU4IN 
/BLCK6/ 
EW,U5IN 


,K34R,K34I,K56R,K56I 

0) ,Y(400) 

(400), U22I(400),V22R<400)  ,  V22K400I  » 

221(400) 

(400) vU33I(400)tV33R(400i ,V33I(400) , 
331(400) 

(400) ,U44I(400),V44R(400) ,  V44 1(400) , 
44I( 400) 

(400) ,U55I(400)  , V5  5R(400)  ,  V55H400)  , 
551(400) 
X,Y,REDFRQ,VRPANL,VIPANL, VRPLNG, V  I PLNG, 

U2  2R,U22I,V2  2R,V2  2I,C22R,C22I, 

EW, V2RNEW, V2 INEW, C2RNEW,C2 INEW 

U33K,U33I,V3  3R, V3  3I,C33R,C33I, 
EW , V3RNE W , V 3  I NE W , C  3RNE  W, C 3 INEW 

IHAVEPtNGRDFNtNUM,NUMBLO,NSTPTStLCOUNTt 
PRES1,IPRES2 

U44R,U44I , V44R, V44I ,C44R,C44I , 
lE  W  ,  V4R  NE  W ,  V4  I  NEW  ,  C4RNE W ,  C  4  1  NE  W 

U55R,U5  5I ,V5  5R, V5  5I ,C55R,C55I, 
EW  ,  V  5R  NE  W , V5 I  NEW  v  C  5RNEW, C SINEW 


J  =  IHAVEP 

I  =  IHAVEP  +  1 

XI  =  X( i) 

ICO   INSURES  PROPER  SELECTION  OF  FLOW  FIELD 
QUANTITIES  AT  F21  IF  THE  GRID  POINT  IS  THE  FIRST  ONE 
AFTER  A  BLADE  OR  WAKE  (  ICO  =  1  ) 

IF  (  ICO  .NE-  1  )  GO  TO  1 

(NEW)  REFERS  TO  THE  PREVIOUS  GRID  POINT  ON  THIS  MACH 
LINE  (IHAVEP).  THIS  GRID  POINT  IS  ON  THE  LOWER  SURFACE 
OF  A  BLADE  OR  WAKE 

A  =  U3RNEW 
B  =  U3INEW 
C  =  C3RNEW 
D  =  C3INEW 
GO  TO  2 

1  A  =  U2RNEW 

B  =  U2INEW 

C  =  C2RNEW 

D  =  C2INEW 

THIS  PORTION  IS  FOR  PITCH 


C      DEFINE  CONSTANTS.  (J)  SUBSCRIPT  KEFERS  TO  PERTURBATION 

C      QUANTITIES  ON  THE  PREVIOUS  MACH  LINE  AT  GRID  POINT 

C      IHAVEP,  (I)  REFERS  TO  QUANTITIES  ON  THE  PREVIOUS  MACH 

C      LINE  AT  (IHAVEP+1),  WHILE  (NEW)  REFERS  TO  QUANTITIES 

C      ON  THE  CURRENT  MACH  LINE  AT  (IHAVEP) 

C 

2  K12R   =    U22RU)     +    C22RU)    4-    AI*U22I(J) 
K12I    =    -AI*U22R(J)     +    U22KJ)    +    C22KJ) 

K34R    =    U22R(I)     -    V22RU)     *■    BI*(     U22KI)     -    C22KI)     ) 

K34I    =    U22KI)    -    V22KI)    -    BI*(     U22R(I)     -    C22R(I)     ) 

K56R    =    A    4-    V2RNEW    4-    BI*(     B    -    D     ) 

K56I    =    B    4-    V2INEW    -    BI*(    A    -    C     ) 
C 

Gl    =    .5*(    K34R    4-    K56R     ) 

G2    =    BI*KL2I 

G3   =    1.    -    AI*BI 

G4    =    ,5*(     K34I    +    K56I     ) 

G5    =    BI*K12R 

G6   =    BI    <-    BI 

G7    =    G3*G3    4-    G6*G6 

G8    =    Gl    -    G2 

G9   =    G4    4-    G5 
C 

C  NOW    STORE    THE    PERTURBATION    QUANTITIES       FOR    GRID    POINT 

C  IHAVEP    ON    THIS    MACH    LINE 

C 

U22RU)     =    U2RNEW 

U22KJ)     =    U2INEW 

V22RU)     =    V2RNEW 

V22HJ)     =    V2INEW 

C22RU)     =    C2RNEW 

C22KJ)     =    C2INEW 
C 

IF    (    ICO    .EQ.    0    )    GO    TO   3 

U33RUJ     =    U3RNEW 

U33HJ)     =    U3IMEW 

V33R(J)     =    V3RNEW 

V33KJ)  =  V3INEW 

C33RU)  =  C3RNEW 

C33KJ)  =  C3INEW 
C 

C      COMPUTE  PERTURBATION  QUANTITIES  AT  THE  CURRENT  GRID 
C      POINT  (IHAVEP4-1),  AND  TEMPORARILY  STORE  AS  (NEW) 

3  U2RNEW  =  (  G8*G3  4-  G9*G6  )/G7 
U2INEW  =  (-G3*G6  4-  G9*G3  )/G7 
V2RNEW  =  .5*(  K56R  -  K34R  ) 
V2INEW  =  .5*1  K56I  -  K34I  ) 
C2RNEW  =  K12R  -  U2RNEW  4-  U2INEW-AI 
C2INEW  =  K12I  -  U2INEW  -  U2RNEW*AI 

C 

C  THIS    PORTION    IS    FOR    PLUNGE.    STEPS    ARE    THE    SAME    AS    FOR 

C  PITCH 


C 


IF    (     ICO    .NE.    1     )    GO    TO    4 

A    =    U5RNEW 

B    =    U5INEW 

C    =    C5RNEW 

D    =    C5INEW 

GO   TO    5 

4  A  =  U4RNEW 
B  =  U4INEW 
C  =  C4RNEW 
D    =    C4INEW 

5  K12R    =    U44R(J)     4-    C44R(J)     4-    AI*U44I(J) 
K12I    =    -AI*J44R(J)     f    U44KJJ     4-    C44KJ) 

K34R  =  U44R(I)     -    V44R(I)     4-  BI*(     1J44K1)     -    C 441(1)     ) 

K34I  =  U44HI)     -    V44iUJ    -  BI*(     U44R(U     -    C44R(I)     i 

K56R  =  A    4-    V4RNEW    4-    BI*(    B  -    D    ) 

K56I  =  B    4-    V4INEW    -    BI*(    A  -    C     ) 


Gl  =  .5*(  K34R 
G2  =  BI*K12I 
G4  =  .5*(  K34I 
G5  =  BI*K12R 
G8  =  Gl  -  G2 
G9  =  G4  «-  G5 


+  K56R  ) 


U44RU) 
U441( J) 
V44R< J) 
V44K  J) 
C44R< J) 
C44KJ) 


U4RNEW 
U4INEW 
V4RNEW 
V4INEW 
C4RNEW 
C4INEW 


IF    (     ICO    .EQ.    0     )    GO    TO    6 

U55RU)  =    U5RNEW 

U55KJ)  =    U5INEW 

V55RU)  =    V5RNEW 

V55KJ)  =    V5INEW 

C55R(J)  =    C5RNEW 

C55KJ)  =    C5INEW 


RESET 
POINT 


THE  LOGIC  VARIABLE  TO  INDICATE  A  GENERAL  FIELD 


ICO  =  0 

U4RNEW  = 

U4INEW  = 

V4RNEW  = 

V4INEW  = 

C4RNEW  = 

C4INEW  = 


(  G8*G3  «■  G9*G6  )/G7 

l-G8*G6  +  G9*G3  >/G7 

.5*(  K5  6R  -  K34R  J 

.5*<  K56I  -  K34J.  ) 

K12R  -  U4RNEW  +  U4INEW*AI 

K12I  -  U4INEW  -  U4RNEW*AI 


RETURN 

FORMATC2H  G,7E16.7,/) 

END 


THIS  IS  SUBROUTINE  NEWBLD 

SUBROUTINE  NEWBLD 

NEWBLD  COMPUTES  U,V  AND  C  AT  THE  FIRST  POSITION  ON  A  6 
ENCOUNTERED  FOR  THE  FIRST  TIME 


REAL*8  X,Y,DSTSTR. 
REAL  K12R,K12IfK34 
DIMENSION  X(400J ,Y 
DIMENSION  U22RC400 

C22R(400).C22 H 
DIMENSION  U33R(400 

C33R(400),C33I( 
DIMENSION  U44R(400 

C44R(400) ,C44I( 
DIMENSION  U55RC400 

C55R<400) ,C55I( 
COMMON  /BLCK1/  X,Y 

S  f  W  f  A I  ,  B I 
COMMON  /BLCK2/  U22 

U2RNEW, U2INEW,V 
CCMMON  /BLCK3/  U33 

U3RNEWf U3IMEWi V 
COMMON  /BLCK4/  I  HA 

JLINEt  ICCIPRES 
COMMON  /BLCK5/  U44 

COM 

USRNEW,U:>li\lEW,  V/ 
COMMON    /GEOM/    FSTR 


HDSTRLtTRNGLHi STAGtFST 

R,K34I,K56R,K56I 

(400) 

i ,U22I(400), V22R(400) , 

400) 

) »U33I(400),V33R(400) , 

400) 

),U44I(400) ,V44R(400) , 

400) 

) ,U55I(400),V55R<400) » 

400) 

»REDFRQ,VRPANLf VIPANLf 

R,U22I,V22R,V22I,C22R, 
2RNEW,V2INEW,C2RNEW,C2 
R,U33I ,V33R»V33I ,C33R, 
~NEW,V3INcW,C3RNEWfC3 
VEPf NGRDFNfNUM,NUMBLDi 
1,1  PRE S2 
R »  U44I  ,  V44R  ,  V44 1  ,  C  44R 

I 

5R  NE  W ,  V  j  1  NEW  .  C  i;RNE  W  ,  C  -j 

MN,DSTSTR,HDSTRL,TRNGL 


RMN 

V22K400)  » 
V33H400)  , 
V44K400)  , 
V55K400)  , 
VRPLNG,VIPLNG, 

C22I, 

INEW 

C33I, 

INEW 

NSTPTS,LCOUNT, 

C44I, 

C55I, 
INEW 
H,STAG,XSUB0, 
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1    DELTA,FAZE 

J  .=  IHAVEP 

I  =  IHAVEP  +  1 

XI  =  X(  I) 

NUM  =  NUM  +  I 

XN  =  FLOAT(NUM) 
C 

C      THIS  PORTION  IS  FOR  PITCH 
C 

C      DEFINE  CONSTANTS.  (J)  SUBSCRIPT  REFERS  TO  PERTURBATION 
C      QUANTITIES  ON  THE  PREVIOUS  MACH  LINE  AT  GRID  POINT 
C      IHAVEP.  (NEW)  INDICATES  QUANTITIES  ON  THE  CURRENT  MACH 
C      LINE  AT  IHAVEP 

K12R  =  U22R(J)  +  C22RU)  +  AI*U22I(J) 

K12I  =  -AI*U22R(J)  +  U22KJ)  +  C22KJ) 

K56R  =  U2RNEW  +  V2RNEW  +  BI*(  U2INEW  -  C2INEW  ) 

K56I  =  U2INEW  +  V2INEW  -  BI*(  U2RNEW  -  C2RNEW  ) 
C 

C      NOW  STORE  THE  PERTURBATION  QUANTITIES   FOR  GRID  POINT 
C      IHAVEP  ON  THIS  MACH  LINE 
C 

U22R(J>  =  U2RNEW 

U22KJ)  =  U2INEW 

V22R(J)  =  V2RNEW 

V22KJ)  =  V2INEW 

C22RU)  =  C2RNEW 

C22KJ)  =  C2INEW 
C 

ARG  =  DELTA*XN 
C 
C      THE  FLOW  TANGENCY  CONDITION 


C 


C 


K34R   =    -(     COS(ARG)     +    XSUBO*REDFRQ*SIN( ARG)     )/S 
K34I    =    -(     SIN(ARG)    -    XSUBO*REDFRQ*COS( ARG)     )/S 


Gl    =    1.    -    AI*BI 

G2   =    BI    *    BI 

G3    =    G1*G1    +•   G2*G2 

G4    =    K56R   -    K34R    -    BI*K12I 

Gi>    =    K56I    -    K34I    «-    BI*K12R 
C 

C  COMPUTE    PERTURBATION    QUANTITIES    AT    THE    CURRENT    GRID 

C  POINT     (IHAVEP+1),     AND    TEMPORARILY    STORE    AS    (NEW) 

C 

U2RNEW    =    (    G4*G1     +    G5*G2    )/G3 

U2INEW    =    (-G4*G2    +    G5*G1    )/G3 

V2RNEW    =    K34R 

V2INEW    =    K34I 

C2RNEW    =    K12R   -    U2RNEW    ♦    AI*U2INEW 

C2INEW    =    K12I    -    U2INEW    -    AI*U2RNEW 
C 

K34R    =    U22R(I)    -    V22R(I)    +    BI*(     U22KI)    -    C22KI)     ) 

K34I    =    U22KI)    -    V22KI)    -    BI*(    U22R(I)    -    C22R(I)     ) 

K56R  =  V2RNEW 

K56I  =  V2INEW 

G4  =  K56R  +  K34R  -  BI*K12I 

G5  =  K56I  +  K34I  *  BI*K12R 

U3RNEW  =  (  G4*G1  +  G5*G2  )/G3 

U3INEW  =  (-G4*G2  +    G5*Gi  ) /G3 

V3RNEW  =  K56R 

V3INEW  =  K56I 

C3RNEW  =  K12R  -  U3RNEW  «■  AI*U3INEW 

C3INEW  =  K12I  -  U3INEW  -  AI*U3RNEW 
C 

C      THIS  PORTION  IS  FOR  PLUNGE.  STEPS  ARE  THE  SAME  AS  FOR 
C      PITCH 


K  1  2P  =  U  44R  (  J  )  *  C  44  R (J )  +  A  I  *  U4  4 1  (  J  ) 

K12I  =  -AI*U44R(J)  +  U44HJ)  *  C44HJ) 

K56R  =  U4RNEW  +  V4RNEW  ♦  BI*(  U4INEW  -  C 41  NEW 

K56I  =  U4INEW  *  V4INEW  -  BI*(  U4RNEW  -  C4RNEW 
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U44RU)     =    U4RNEW 

U44KJ)     =    U4INEW 

V44RU)     ■    V4RNEW 

V44KJ)     =    V4INEW 

C44RU)     =    C4RNEW 

C44KJ)     =    C4INEW 

K34R    =    REDFRQ*SIN( ARG)/S 

K34I    =    -    REDFRQ*COS(ARG)/S 

G4    =    K56R    -    K34R    -     BI*K12I 

G5    =    K56I    -    K34I     f    BI*K12R 

U4RNEW    =    (     G4*G1     +    G5*G2    )/G3 

U4INEW    =    (-G4*G2    ♦    G5*G1     )/G3 

V4RNEW    =    K34R 

V4INEW    =    K34I 

C4RNEW    =    K12R   -    U4RNEW    +    AI*U4INEW 

C4INEW    =    K12I    -    U4INEW    -    AI*U4RNEW 

K34R    =    U44R(I)     -    V44R(I)     +    BI*(     U44KI) 

-    C44K  I  )     ) 

K34I    =    U44HI)     -    V44KI)    -    BI*(     U44RU) 

-    C44R(I)     ) 

K56R   =    V4RNEW 

K56I    =    V4INEW 

G4    =    K56R    +    K34R    -    BI*K12I 

G5    =    K56I     *    K34I     +    BI*K12R 

U5RNEW    =    I    G4*G1    +    G5*G2    )/G3 

U5INEW    =    <-G4*G2    +    G5*G1    )/G3 

V5RMEW    =    K56R 

V5INEW    =    K56I 

C5RNEW    =    K12R   -    U5RNEW    *    AI*U5INEW 

C5INEW    =    K12I    -    U5INEW    -    AI*U5RNEW 

NNBLD    =    NUM    +    1 

RETURN 

00    FORMAT ( 2H    N,7E16.7,/) 

01    FORMAT < s0« ,10X,8H    BLAOE    #,13) 

END 

THIS  IS  SUBROUTINE  WAKE 

SUBROUTINE  WAKE 

WAKE  COMPUTES  U,V  AND  C  ABOVE  AND  BELOW  THE  SLIP  PLANE 
AFT  OF  A  BLADE. 

REAL*8  X,Y 

DIMENSION  X(40 

DIMENSION  U22R 
1    C22R(400),C 

DIMENSION  U33R 
I    C33R(400),C 

DIMENSION  U44R 
1    C44R(400),C 

DIMENSION  U5  5R 
1    C55R(400),C 

DIMENSION  A(8, 

COMMON  /BLCK1/ 
1    S,W,AI,BI 

COMMON  /BLCK2/ 
1    U2RNEWfU2IN 

COMMON  /BLCK3/ 
1    U3RNEW,U3IN 

COMMON  /BLCK4/ 
i    JLINE,ICO,I 

COMMON  /BLCK5/ 
1    U4RNEW.U4IN 

COM'1  ON  /BLCK6/ 
1    U5RNEW,U5IN 

J  =  IHAVEP 


D0)f V22R(400) ,V22I(400) , 
30),V33R(400)  ,  V33K400)  , 
)0),V44R<400) , V44I1400) , 
30),V5  5R(400) ,V 55  I  (400), 


0)  ,Y(400) 

(400) ,U22I(400 

221(400) 

(400)  .U33K40C 

331(400) 

(400)  ,U44I(40C 

441(400) 

(400) fU55I(40C 

55  1(400) 

x!y,redfrq,vrpanl,vipanl,vrplng,viplng, 

u22r,u2  2i ,v2  2r, v2 21, c22r,c 22  i, 

e  w  ,  v  2r ne h ,  v2  i  ne  w  ,  c  2rne  w ,  c  2  i  ne  w 

U33R,U33I,V3  3R,V3  3I,C33R,C33If 

EW,V3RNEW»V3lNEW,C3RNEW,C3INEW 
I H AVER, NGRDFN,NUM,NUMBLD,NSTPTS,L COUNT, 

PRESi, IPRES2 
U44R,U44I,V44R,V44I,C44R,C44I, 

EW,V4RNEW,V4INEWtC4RNEW,C4INEW 
U55R.U55I ,V55P , V5  5 ;  ^i, 

EWtV5RNEWtV5INtW»C5RNEW,C        NEW 
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I  =  IHAVEP  +  1 
XI  =  X(  II 

ITIME  TELLS  WHETHER  PLUNGE  SOLUTION  HAS  BEEN  DONE.  IF 
ITIME  =  2,  PLUNGE  SOLUTION  IS  COMPLETE;  IF  ITIME  =  1, 
PLUNGE  SOLUTION  NEEDS  TO  BE  DONE 

ITIME  =  1 

THIS  PORTION  IS  FOR  PITCH 

COMPUTE  THE  RIGHT  HAND  SIDE  OF  THE  EQUATIONS.  (I)  AND 
(J)  SUBSCRIPTS  REFER  TO  PERTURBATION  QUANTITIES  AT 
GRID  POINTS  (IHAVEP+1)  AND  (IHAVEP)  ON  THE  PREVIOUS 
MACH  LINE,  WHILE  (NEW)  REFERS  TO  GRID  POINT  (IHAVEP) 
ON  THE  CURRENT  MACH  LINE 


B(l) 
B(2) 
B(3) 
B(4) 
B(5) 
6(6) 
B(7) 
B(8) 

NOW  ST 
IHAVEP 

U22R( J) 
U22K  J) 
V22R( J) 
V22KJ) 
C22R( J) 
C22KJ) 


U22RU)    +    C22R(J) 
U22KJ)    *-    C22KJ) 
U2RNEW    +•    V2RNEW 
U2INEW    *    V2INEW 


U33R( J) 
U33K  J) 
U22R(  I) 
U22K  I) 


4-  C33R(J) 


*■    AI*U22I(J) 
-  AI*U22R(J) 
3I*(  U2INEW  - 
BI*(  U2RNEW  - 
AI-U33KJ) 


C2INEW  ) 
C2RNEW  ) 


+  C33KJ)  -AI*U33R(J) 


V22R(  I) 
-  V22KI) 


BI*( 

BI#( 


U22K  I) 
U22RU) 


-  C22K  I  ) 

-  C22R(I) 


ORE  THE  PERTURBATION  QUANTITIES 
ON  THIS  MACH  LINE 

U2RNEW 
U2INEW 
V2RNEW 
V2INEW 
C2RNEW 
C2INEW 


FOR  GRID  POINT 


SET  UP  THE  MATRIX  OF  COEFFICIENTS 

DO  2  K 

DO  2  L  =  1,8 

A(K,L) 

A(l,l) 

A(l,7) 

A(2,2) 

A(2,8) 

A(3,l) 

A(3,5) 

A(4,2) 

A(4,6) 

A(5,3) 

A(5,7) 

A(6,4) 

A(6,8) 

A (7,3) 

A(8,4) 

A(l,2)  =  -AI 

A(5,4)  =  -AI 

A(6,3) 

A(2,l) 

A(3,2)  =  -BI 

A(4,7)  =  -BI 

A(7,4)  =  -BI 

A(8,7)  =  -BI 

A(3,8)  =  BI 

A(4,l) 

A(7,8) 

A  (  8 ,  3  ) 

A(7,5) 

A(8,6) 

CALL  SIMULTANEOUS  LINEAR  EQUATION  SOLVER  TO  CALCULATE 
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c 
c 

c 

U,  V  AND  C   AT  THIS  GRID  POINT  (IHAVEPH) 

CALL  SIMQ< A,B,8,KS) 

IF  (  ITIME  .EQ.  2  )  GO  TO  3 

c 

ITIME  =  2 

U2RNEW  =  B(l) 

U2INEW  =  B(2) 

U3RNEW  =  B(3) 

U3INEW  =  B(4) 

V2RNEW  =  B<5) 

V2INEW  =  6(6) 

C2RNEW  =  B(7) 

C2INEW  =  B(8) 

V3RNEW  =  B(5) 

V3INEW  =  8(6) 

C3RNEW  =  B(7) 

C3INEW  =  B(8) 

c 
c 

COMPUTE  THE  RIGHT  HAND  SIDE  OF  THE  EQUATIONS 

c 
c 

THIS  PORTIGN  IS  FOR  PLUNGE.  STEPS  ARE  THE  SAME  AS  F 

c 

c 

PITCH 

B(l)  =  U44RU)  +  C44R(J)  *-  AI*U44I(J) 

B(2)  =  U44HJ)  ♦  C44KJ)  -  AI*U44R(J) 

B(3)  =  U4RNEW  *•  V4RNEW  +  BI*(  U4INEW  -  C4INEW  ) 

B(4)  =  U4INEW  *  V4INEW  -  BI*(  U4RNEW  -  C4RNEW  ) 

B(5)  =  U55RU)  *  C55RU)  +  AI*U55I(J) 

B(6)  =  U55HJ)  +  C55KJ)  -AI*U55R(J) 

8(7)  =  U44P(I)  -  V44R(I)  ♦  BI*(  U44IU)  -  C 441(1)  ) 

c 

B(8)  -  U44KI)  -  V44KIJ  -  BI*(  U44R(I)  -  C44R(  I  )  ) 

U44R(J)  =  U4RNEW 

U44KJ)  =  U4INEW 

V44RU)  =  V4RNEW 

V44HJ)  =  V4INEW 

C44R(J)  =  C4RNEW 

c 
c 
c 

c 

C44KJ)  =  C4INEW 

SET  UP  THE  MATRIX  OF  COEFFICIENTS 

GO  TO  1 

3  U4RNEW  =  B(l) 

U4INEW  =  B(2  J 

U5RNEW  =  B(3) 

U5INEW  =  B(4J 

V4RNEW  =  B(5) 

V4INEW  =  B(6) 

C4KNEW  =  B(7) 

C4INEW  =  B(8) 

V5RNEW  =  B(5) 

V5INEW  =  B(6) 

C5RNEW  =  B(7) 

C5INEW  =  B(8) 

ITIME  =  1 

RETURN 
100  F0RMATC2H  W.7E16.7,/) 

END 


C  THIS  IS  SUBROUTINE  SIMQ 

C      THIS  IS  SUBROUTINE  SIMQ 

SUBROUTINE  SIMQ( A, b» N, KS) 
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C  SIMQ  IS  A  SIMULTANEOUS  LINEAR  EQUATION  SOLVER 

C 

C  A   -   MATRIX  (N  BY  N)  OF  COEFFICIENTS  STORED 

C  COLUMNWISE  (DESTROYED) 

C  B   -   VECTOR  OF  ORIGINAL  CONSTANTS.   FINAL  SOLUTION 

C  ON  THIS  VECTOR. 

C  N   -   NUMBER  OF  EQUATIONS  AND  VARIABLES 

C  MUST  BE  >  1 

C  KS  -   OUTPUT  DIGIT 

C  0   FOR  NORMAL  SOLUTION 

C  J.   FOR  A  SINGULAR  SOLUTION 


DIMENSION  Al 1) ,8(1) 
C 

C      FORWARD  SOLUTION 
C 

TOL  =  0. 

KS  =  0 

JJ  =  -N 

DO  65  J  =  1,N 

JY  =  J  *  1 

JJ  =  JJ  ♦  N  +  1 

BIGA  =  0. 

IT  =  JJ  -  J 

DO  30  I  =  J, N 
C 

C      SEARDH  FOR  MAXIMUM  COEFFICIENT  IN  COLUMN 
C 

'   IJ  =  IT  +  I 

IF  (ABS(BIGA)  -  ABS(AUJ))  )  20,30,30 
20  BIGA  =  A( IJ) 

I  MAX  =  I 
30  CONTINUE 

C 

C      TEST  FOR  PIVOT  LESS  THAN  TOLERANCE  (SINGULAR  MATRIX) 

C 

IF  (ABS(B-IGA)  -  TOL)  35,35,40 
35  KS  =  1 

RETURN 
C 

C      INTERCHANGE  ROWS  IF  NECESSARY 
C 

40  II  -  J  *  N*(J-2) 

IT  =  I  MAX  -  J 

DO  50  K  =  J,N 

II  =  II  4-  N 
12  =  II  +  IT 
SAVE  =  AII1) 
A(I1)  =  AC  12) 
A(I2)  =  SAVE 

C 

C      DIVIDE  EQUATION  BY  LEADING  COEFFICIENT 

C 

50  A(I1)  -    A( II )/BIGA 

SAVE  =  B( I  MAX) 

B(IMAX)  =  B(J) 

B(J)  =  SAVE/BIGA 
C 

C      ELIMINATE  THE  NEXT  VARIABLE 
C 

IF  (J-N)  55,70,55 
55  IQS  =  N*(J-I) 

DO  65  IX  =  JY,N 

IXJ  =  IQS  ♦  IX 

IT  =  J  -  IX 

DO  60  JX  =  JY,N 

IXJX  =  N*(JX-1)  ♦  IX 

JJX  =  IXJX  f  IT 
60  A(IXJX)  =  A(IXJX)  -  (A(IXJ)*A( JJX) ) 
6  5  B(IX)  =  BUX)  -  (BU)*A(  IXJ)  ) 
C 
C      BACK  SOLUTION 
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70  NY  =  N-l 

IT  =  N*N 

DO  80  J  =  1,  NY 

IA  =  IT  -  J 

IB  =  N  -  J 

IC  =  N 

DO  30  K  =  I, J 

B(I6)  =  B( IB)  - 

A(  IA)*B( IC) 

IA  =  IA  -  N 

80  IC  =  IC  -  1 

RETURN 

END 

THIS  IS  SUBROUTINE  COEF 

SUBROUTINE  COEF 

SUBROUTINE  COEF  COMPUTES  THE  NON-DIMENSIONAL  LIFT 
FORCE  AND  MOMENT  FOR  EACH  BLADE.  IN  ADDITION,  THE 
MAGNITUDE  AND  PHASE  OF  THE  COMPLEX  LIFT  FORCE  AND 
MOMENT  ARE  DETERMINED 


REAL* 8 
1    XN 

DIMENSI 
1  100) 
2100) ,RM 

DIMENSI 

1  RM2( 

2  ANGP 

3  R12M 
DIMENSI 
DIMENSI 

1  C22R 
DIMENSI 

1  C33K 
DIMENSI 

1  C44R 
DIMENSI 

1  C55R 
DIMENSI 
COMMON 

1   s,w, 

COMMON 
1    U2RN 

COMMON 
1    U3RN 

COMMON 
1    JLIN 

COMMON 
1    U4RN 

COMMON 
1    U5RN 

COMMON/ 
1    ANGM 

COMMON/ 
1    NITO 

COMMON 
1    DELT 

DATA  LB 

1  4,1 

2  1,2 

3  2,3 


X,Y,FSTRMN,DSTSTR,HDSTRL,TRNGLH,STAG,XNEW,XOLD, 


ON    XPT(4, 100)  ,PR12(4,  100) ,PIi2(4 

,RM 112(4, 100) ,PR34(4,100)  ,  PI  34(4 

134(4,100) 

ON    RLK  76)  ,RL2(76)  ,RL3(  76  ),RL4(7 

76) ,RM3( 76) ,RM4( 76) , ANGP12C 76 ) , A 

34( 76) , ANGM34( 76) ,P12MAG( 76) ,P34 

AG( 76),R3  4MAG(76) 

ON    X(400) ,Y(400) 

ON    U22R(400i ,U22I(400) , V2  2R(400) 

(400) ,C22I(400) 

ON  U33R(4  00) ,U33I ( 400) , V33R(400 ) 

(400) ,C33I (400) 

ON  U44R(400) ,U44I(400) ,V44R(400) 

(400) ,G44I(400) 

ON  U5  5R(400) , U5  5 I ( 400 ) , V55R(400) 

(400) ,C55I(400) 

ON  NN(IOO) ,10(76) ,LB(76) 

/BLCK1/  X,Y,REDFRQ,VRPANL,VIPANL 

AI  ,BI 

/BLCK2/  U22R,U2  2I ,V2  2R, V2  2I,C22R 

EW, U2INEW,V2RNEW, V2INEW,C2RNEW,C 

/BLCK3/  U33R,U3  3I,V3  3R, V3  3I,C3  3R 

EW,U3INEW,V3RNEW,V3INEW,C3RNEW,C 

/BLCK4/  IHAVEP,NGRDFN,NUM,NUMBLD 

E,IC0,IPRES1,IPRES2 

/BLCK5/  U44RTU44I,V44R,V44I,C44R 

EW,U4INEW,V4RNEW,V4INEW,C4RNEW,C 

/BLCK6/  U55R,U55I,V5  5R, V5  5I,C55R 

E W , U5 I NE W , V  5R  NE  W , V5 I NE  W , C  5RNE  W , C 

BLCK7/  RL1 , RL2, RL3, RL4, RM1,RM2,R 

12, ANGP34, ANGM34,P12MAG,R12MAG,P 

BLCKO/    XPT,PR12,PU2,PR34,PI34,  I 

T,NJUNK 

/GEOM/    FSTRMN,DSTSTR,HDSTRL,TRNG 

A, FAZE 

/It  2, 3, 4,1, 2, 3, 4, 1,2, 3, 4,  1,2, 3, 4 

2,3,4,1,2,3,4,1,2,3.4,1,2,3,4,1, 

3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2, 

4/ 


,100) ,RMR12(4, 
,  100) ,RMR34(4, 

6),RM1(76) , 
NGM12( 76)  , 
MAG( 76), 


,V22I(400)  , 
,  V33K400)  , 
,  V44K400)  , 
,V55I(400) , 

,VRPLNG,VIPLNG, 

,C22I, 

2  I  NEW 
, C33I, 

3INEW 
,NSTPTS,LCCUNT, 

,C44I, 

4INEW 

,C55I, 

5  I  MEW 

M3,RM4,ANGP12, 

34MAG,R34MAG 

WRITE, LVEC, 

LH,STAG,XSUBO, 

,1,2,3,4,1,2,3, 

2,3,4,1,2,3,4, 

3,4,1,2,3,4,1, 


ENDPT(Y2,Y1,X0,X2,DELX) 

1         X2    ) /DELX    +    Y2 


(    Y2    -   Yl    )*(    XO    - 
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XARM(V)  =  V  -  XSUBO 
C 
C 

C  INITIALIZE  THE  ARRAYS  AND  SET  VALUES  AT  THE  LEADING 

C  EDGE  OF  THE  FIRST  BLADE.  CALLED  AFTER   INITAL.   ID(L) 

C  IS  A  BLADE  DEPENDENT  SWITCH  VARIABLE: 

C         =0   GRID  POINT  IS  ON  THE  BLADE  —  MAKE  CALCULATION 
C         =1   GRID  POINT  IS  IN  THE  WAKE  —  SKIP  IT 


C 


C 


DO  10  IL  =  1,4 
DO  5  IJ  =  It  100 
NN(IJ)  =  0 
XPT( IL, IJ)  =  0.0 
PR12UL,  IJ)  =  0. 
PI12(IL,I J)  =  0. 
PR34UL,  IJ)  =  0. 
PI34<  ILtU)  =  0. 
RMR12( I L, IJ)  =  0. 
RMI12( ILt IJ)  =  0. 
RMR34( IL, I  J)  =  0. 
5  RMI34( IL, I  J)  =  0. 

10  CONTINUE 

DO  11  IL  =  1,76 

11  ID(IL)  =  0 


NN(1)  =  1 

Al  =  0.5*(  C55R(1)  -  C44R(1)  ) 

A2  =  0.5*(  C55K1)  -  C44K1)  ) 

A3  =  0.5*1  C33R(1)  -  C22RI1)  ) 

A4   =    0.5*(    C33H1)    -    C22U1)     ) 

PR  12 (1,1)  =  Al 

PI12(1,1)  =  A2 

PR34(1,1)  =  A3 

PI34( 1,1)  =  A4 

RMR12C 1,1)  =  -  XSUB0*A1 

RMI12(1 ,1)  =  -  XSUB0*A2 

RMR34(1,1J  =  -XSUB0*A3 

RMI34( 1,1)  =  -XSUB0*A4 

XPT(1,1)  =  SNGL(  X(l)  ) 

RETURN 
C 

ENTRY  C0EF1 
C 

C      THIS  SECTION  CALCULATES  DELTA  PRESSURE  AND  DELTA 
C      MOMENT  FOR  EACH  GRID  POINT  ON  THE  BLADE.  CALLED  AFTER 
C      SUBROUTINES  BOTTOM  OR  NEWBLD. 
C 

M  =  NUM  +  1 

I  =  IHAVEP  +  1 

L  =  LB(M) 

NN(L)  =  NN(L)  ♦  1 

NI  =  NN(L) 

XN  =  DFLOATI  NUM  ) 

XNEW  =  X( I )  -  STAG*XN 

XPT(L,NI)  =  SNGL(  XNEW  ) 

Al  =  0.5*(  C5RNEW  -  C4RNEW  ) 

A2  =  0.5*(  C5INEW  -  C4INEW  ) 

A3  =  0.5*(  C3RNEW  -  C2RNEW  ) 

A3  =  0.5*(  C3INEW  -  C2INEW  ) 
C 

C      ARG1  IS  USED  TO  COMPUTE  THE  PRESSURE  ON  THE  N(TH) 
C      BLADE  WHEN  IN  THE  INITIAL  POSITION. 
C 

C  ARG2  IS  USED  TO  COMPUTE  THE  PRESSURE  ON  THE  N(TH) 
C  BLADE  WHEN  IN  THE  MEAN  POSITION,  PITCHING  LEADING 
C      EDGE  UP. 


ARG1  =  -  DELTA*SNGL(XN) 

ARG2  =  4.712389  -  DELTA*SNGL(XN) 

CARG1  =  COS(ARGl) 

CARG2  =  C0SIARG2) 
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SARG1    =    SIN(ARGl) 

SARG2    =    SINCARG2) 
C 

PR12<L,NI)    =    A1*CARG1    -    A2*SARG1 

PI12(LfNII    =    A1*CARG2    -    A2*SARG2 

PR34<L,NI)     =    A3*CARG1    -    A4*SARGl 

PI34(LtNII    =    A3*CARG2    -    A4*SARG2 

RMR12CL.NI)    =    PR12(LtNI)*XARM(XNEW) 

RMU2(L,NI)    =    PI12(L,NI)*XAkM(XNEW) 

RMR34(L,NI)  =  PR34 < L, N I ) *XARM (XNEW J 

RMI34(L,NI)  =  PI34(L, NI ) *XARM< XNEW) 

GO  TO  35 
C 

ENTRY  COEF2 
C 

C      FOR  A  PARTICULAR  BLADE.   THE  VALUES  OF  DELTA  PRESSURE 
C      AND  MOMENT  ARE  EXTRAPOLATED  LINEARLY.  THIS  SECTION  IS 
C      CALLED  AND  USED  THE  FIRST  TIME   WAKE   IS  USED  FOR  EACH 
C      BLADE. 
C 

M  =  NUM  +  1 

IWRITE  =  0 
C 

C      IF  THIS  IS  NOT  THE  FIRST  CALL  FROM  THE  WAKE  AFT  OF 
C      THIS  BLADE,  RETURN 
C 

IF  (  ID(M)  .EQ.  1  )   GO  TO  35 

I  =  IHAVEP  +  1 

L  =  LB(M) 

NN(L)  =  NN(L)  +  1 

NI    =  NN(L) 

NIM  =  NI  -  1 

ID(M)  =  1 
C 

C      THESE  THREE  VARIABLES  ARE  USED  IN  THE  MAIN  PROGRAM  IF 
C      THE  PRESSURE  DISTRIBUTION  IS  OUTPUT. 

C         IWRITE  =  1  —  CALCULATIONS  FOR  BLADE  ARE  COMPLETE 
C         NITOT   —   TOTAL  NUMBER  OF  GRID  POINTS  ON  BLADE 
C         LVEC   —   ENTRY  POINT  FOR  THE  STORAGE  ARRAYS 


C 


C 


IWRITE  =  1 
NITOT  =  NI 
LVEC    =    L 


XN    =    DFLOAT(NUM) 

XNEW  =  X<I  )  -  STAG*XN 

XPT(L,NI)  =  1.0 

XPT(L,1 )  =  0. 

XOLD  =  XNEW  -  DSTSTR 
C 
C      TRAPEZOIDAL  INTEGRATION  FOR  LIFT  FORCE  AND  MOMENT 


SUM1  =  0. 

SUM2  =  0. 

SUM3  =  0. 

SUM4  =  0. 

SUM5  =  0. 

SUM6  =  0. 

SUM7  =  0. 

SUM8  =  0. 
NI2  =  NI  -  2 
DO  71  IJ  =  2,NI2 

SUM1  =  SUM1  +  PR12(LfIJ) 

SUM2  =  SUM 2  *  PI  12 (L, IJ  1 

SUM3  =  SUM3  ♦  RMR12(L,IJ) 

SUM4  =  SUM4  +    RMI12(LtIJJ 

SUM5  =  SUM5  +  PR34(L,IJ) 

SUM6  =  SUM6  +  PI34(L,IJ) 

SUM7  =  SUM7  +  RMR34(L,IJ) 

71  SUM8  =  SUM8  +  RMI34(L,IJ) 


S  U  M 1  +  0«5*(  PR12(L,1)  +  PR12(L,NIM)  »  »* 
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SUM6 
SUM7  =  ( 


OSTSTR 
SUM2  =  (  SUM2  +  0.5*1  PI12(L,1)  ♦  PI121L,NIM)  )  )* 

DSTSTR 
SUM3  =  1  SUM3  «■  0.5*1  RMR121L»1)  +  RMR12lL,NIM)  )  )* 

DSTSTR 
SUM4  =  (  SUM4  +  0.5*(  RMIi2(L,l)  +  RMIL21L.NIM)  )  )* 

DSTSTR 
SUM5  =  (  SUM5  *■    0.5*1  PK3<t(L,l)  +  PR341L,NIM)  )  )* 

DSTSTR 
(  SUM6  +  0.5*1  PI34iL,l)  <-  PI34lL,NIM)  )  )* 

DSTSTR 
SUM7  ♦  0.5*1  RMR34lL,l)  *  RMR341L,NIM)  )  )* 

DSTSTR 
SUM8  =  (  SUM8  *■    0.5*(  RMI341L.1)  ♦  RMI34(LtNIMJ  )  )* 

DSTSTR 

EXTRAPOLATE  FOR  VALUES  AT  THE  END  POINT  OF  THE  BLADE. 

PR121L,NI)    =    ENDPT(PR12(L,NIM) tPR12(L,NI2)t l.tXOLO, 

DSTSTR) 
PI12(L,NI)    =    ENDPT(PU2(LtNIM)  ,PI12(LtNI2),l.,XOLOt 

DSTSTR) 
PR34(L»NI)  =  ENDPT1PR34<L,NIM) ,PR34<L,NI2)fl.,X0LDi 

DSTSTR) 
PI34(L,NU  =  ENDPT(PI34(LfNIM) iPI34(LfNI2)i 1-fXOLDf 

DSTSTR) 

RMR12(LfNII    =    PR12(L,NI)*< 1.  -    XSUBO) 

RMI12(L,NI)    =    PU2(L,NI)*<  1.  -    XSUBO) 

RMR34(LfNI)    =    PR341 L, NI )*1 1 .  -    XSUBO) 

RMI341L,NI)    =    PI34(L,NI)*(1.  -    XSUBO) 
C    =    (     1.0    -    XOLD     »*0.5 

RLKM)  =  SUM1  4-    (  PR12(L,NIJ    +    PR12(L,NIM)     )*C 

RL2(M)  =  SUM2  +    1  PI12(L,NI)    +    P 1 12(  L»  N I M)     )*C 

RMUM)  =  SUM3  ♦    (  RMR12(L,NI)     +    RMR121L,NIM)     )*C 

RM21M)  =  SUM4  +    (  RMI12(LtNI)    +    RMI121L,NIM)     )  *C 

RL3(M)  =  SUM5  *     (  PR341L,NI)     *     PK34(LfNlM)     )*C 

RL4(M)  =  SUM6  «■    (  PI34(L,NI*    +    PI34lLfNIM>     )*C 

RM31M)  =  SUM7  +     (  RMR34(L,NI)     +    P.MR341  L  t  NIM  )     )*C 

RM41M)  =  SUM8  *     (  RMI341L»NI)     *    RMI34(L,NIM)     )*C 

PHASE    ANGLE 

ANGP121M) 
ANGM121M) 
ANGP341M) 
ANGM341M) 

MAGNITUDE 

P12MAG1M)  =  SQRT(  RLl {M)*RL1( Mi  +  RL21M) *RL2(M)  ) 

R12MAG1MI  =  SQRT1  RM1 ( M) *RM1 (M )  +  RM21 M) *RM21 M)  ) 

P34MAG1M)  =  SQRT(  RL3( M)*RL3< M J  *  RL4CM) *RL4(M)  ) 

R34MAG1M1  =  SQRT(  RM3( M)*RM3(M )  +  RM4 ( MJ *RM4l M)  ) 


(  ATAN2(  RL2(M),RL1(M)  )  )*57*29578 

(  ATAN21  RM21M) ,RM1(M)  )  )*57. 29578 

(  ATAN2(  RL4(M) tRL3(M)  )  )*57. 29578 

(  ATAN21  RM4(M) ,RM3(M)  )  )*57. 29578 


NN(L) 


35  RETURN 
END 


0 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C  CYLINDRICAL  SHELL/PANEL  FLUTTER  PROGRAM  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  THIS  IS  THE  MAIN  PROGRAM 


C 


COMPLEX    PHI  (400}  .DPHISH400)  ,DPHIS2<400) 

DIMENSION    X(400)  ,R(400)  ,  XT' 1(400)  ,CPR<400)  ,CPI(400) 

DIMENSION  DUMMY( 400) ,DUNCE( 400) 

C0MM0N/BLCK1/  FSTRMN, DELTAS , DX, DH, DSTSTR ,R0  ,RI ,REDFRQ 

COMMON/ BLCK2 /  A  1 , A  2 , A3  T  A4 , RN, RM , RM 1 , F AC  T 1 , F  ACT  2 

COMMON/ BLCK3/  NGRDFN, NPTS , LCOUNT , I HA VEP , N, M ,M1 , 1  PR  INT 

C0MM0N/BLCK4/  PHI , DPH I  SI » DPHI S2 

C0MM0N/BLCK5/  X» R T XPT T CPR , CPI 

OATA  PI/3.141593/ 

ENDPT(Y2,Yi,X0,X2,DELX)  =  (  Y2  -  Yl  >*(  XO  -X2  ) /DELX 
1    +  Y2 


EPS  =  l.E-04 
IER  =  0 

CALL  INPUT(IER) 
IF  (  IER  .EQ.  1  )  GO  TO  12 
1000  CONTINUE 
C 

IHAVEP  =  1 
LCOUNT  =  1 


CALL    INITAL 
1    CONTINUE 
CALL    COMPXR 


C 

C      IS  THE  GRID  POINT  ON  THE  INITIAL  RIGHT  RUNNING  MACH 

C      LINE? 

C 

IF  ((  LCOUNT  .EQ.  1  )  .AND.  (IHAVEP  .LE.  NGRDFN  )) 
1    GO  TO  3 
C 

C      IS  THE  GRID  POINT  ON  THE  OUTER  CYLINDER  SURFACE? 
C 

IF  (  IHAVEP  .EQ.  0  )  GO  TO  5 

C      IS  THE  GRID  POINT  ON  THE  INNER  CYLINDER  SURFACE? 
C 

IF  (  IHAVEP  .EQ.  NGRDFN  )  GO  TO  6 
C 

C      THEN  THE  GRID  POINT  IS  IN  THE  GENERAL  FLOW  FIELD. 
C 

2  CALL  GENFPT( IER) 

IF  (  IER  .EQ.  1  )  GO  TO  12 
IHAVEP  =  IHAVEP  *■    1 
GO  TO  1 
C 

3  CALL  MACHLN 

IF  (  IHAVEP  .EQ.  NGRDFN  )  GO  TO  4 

IHAVEP  =  IHAVEP  +  1 

GO  TO  1 
C 
C      THIS  IS  THE  LAST  POINT  ON  THE  INITIAL  MACH  LINE. 


C 


4  LCOUNT  =  LCOUNT  +  1 
IHAVEP  =  0 
GO  TO  1 


5  CALL  RAD2 
C 
C      IS  THIS  GRID  POINT  PAST  THE  FLEXIBLE  CYLINDER  LENGTH? 
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IF  {  (XI1UEPSI  .GE.  1.  *  GO  TO  8 

IHAVEP  =  IHAVEP  «-  1 
GO  TO  1 

6  CALL  RAD1 
IHAVEP  =  0 
LCOUNT  =  LCOUNT  +  1 
GO  TO  1 

8  CONTINUE 

IF    (     AbS(XPT(LCOUNT)     -    1.     )     .LE.     EPS    )     GO    TO    10 

L    =    LCOUNT 

LM1    =    LCOUNT    -    1 

LM2    =    LCOUNT    -    2 

XPT(L)     =    1. 

CPR(L)     =    ENDPT(CPR(LM1),CPR(LM2)  ,  1.00,  XPHLMi  ), DSTSTR) 

CPKL)     =    ENOPTCCPI  (LM1)  ,CPI(LM2)  f  1.00,  XPHLMi)  ,  DSTSTR) 

WRITE<6,105) 

DO    9    I    =    1,L 

ARG    =    RM1*PI#XPT(I  ) 

WRITE (6, 110)     I,XPT(I),CPR(I),CPI(I) 

DUMMY(I)    =    CPR( I)*SIN( ARG) 

9  OUNCE(I)    =    CPI ( I)*SIN(ARG) 

SUM1    =    0.5=MCPR(L)     +    CPR(LMl) )*( 1.00    -    XPT(LMD) 
SUM2    =    0.5*(CPI(L)     *    CPKLMU  )*(  1.00    -    XPT(LMD) 

WRITE (6, 114)    FSTRMN,REDFRQ,RO,RI,N 

C  ALL    QSF( DSTSTR, DUMMY , DUMMY , LM1 } 
CALL    QSF( DSTSTR, DUNCE, DUNCE, LM1) 
QREAL    =    0.5*(DUMMY(LM1)     ♦    SUMi) 
Q1MAG    =    0.5*(DUNCE(LM1)     «-    SUM2) 

WRITE (6, 11 5)     M,M1,QREAL,M,M1,QIMAG 
GO   TO    12 

10  CONTINUE 
XPT(L)     =    1. 
WRITE(6,105) 

DO    11    I    =    1,  LCOUNT 

ARG    =    RM1*PI*XPT<  I  ) 

WRITE (5, 110)     lyXPT(I) ,CPR(IJ ,CPI(I) 

DUMMY(I)     =    CPR( I )*SIN( ARG) 

11  DUNCE(I)     =    CPIU)*SIN(ARG) 


WRITE (6, 11 4)    FSTRMN, REDFRQ,RO,R I ,N 

CALL    QSF( DSTSTR, DUMMY, DUMMY, L) 

CALL    QSF(DSTSTR, DUNCE, DUNCE, LJ 

QREAL    =    0.5*DUMMY(L) 

QIMAG=    0.5*DUNCE(L) 

WRITE (6, 11 5)  M,M1, QREAL, M, Ml , QI MAG 

12  CONTINUE 
STOP 
150  F0RMAT(F10.5) 

105  FORMAT ( « 1 • , Tl 9, • I • , T2  7 , • X ' , T38 , • CPR* , T5  0 , • C PI • , // ) 
110  FORMAT ( «0«,15X,I5,3(3X,F9.5)) 

114  FORMAT(//«  •,T20,,MACH  NUMBER  =  ',F10.7,//'  '  ,T20, 

1  'REDUCED  FREQUENCY  =  '^10,7,/'  ■ ,T20, M RADIUS/ ■ , 

2  'LENGTH)  OUTER  =  ' ,F10.7,/«  «,  T20 ,«{ RADI US/LENGTH) • 

3  ,«  INNER  =  ',F10.7,/«  ', T20, "CIRCUMFERENTIAL  MOOES 

4  /•  ' ,T34, 'NUMBER  =  ',15,//) 

115  FORMATS '0' ,T20,'Q( • ,12,' ,', 12,' )REAL  =  • ,F10. 5, / • 0 • , 
1    T20,'Q(  ',12,'  ,'  ,I2,M  IMAG  =  «,F10.5,/J 

136  FORMAT! •  !•  ) 
END 
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C  THIS  IS  SUBROUTINE  INPUT 

C 

SUBROUTINE  INPUTUER) 
C 

C      SUBROUTINE  INPUT  READS  INPUT  VARIABLES,  AND  DEFINES 
C      SEVERAL  CONSTANTS  USED  THROUGHOUT  THE  PROGRAM. 


C 


C 


COMMON/ BLC.K1/  FSTRMN,  DELTA  S  ,  DX,  DM  ,  DSTSTR  ,  RO  ,R  I  ,  REDFRQ 
COMMON/BLCK2/  A 1, A2, A3.A4 , RN , RM, RM1. FACT  1 f F ACT2 

COMMON/BLCK3/  NGRDFN , NPTS, LCOUNT , I  HA VEP ,NtM , Ml f I  PRINT 
NAMELIST/NAM1/  FSTRMN, REDFRQ ,RO ,R I , NGRDFN , N , M, Ml , 
1    IPRINT 

WRITE(6,80) 

READ(5, NAM1) 
WRITE(6,NAMi) 

FMANGL    =    ARSINd  ./FSTRMN) 
DH    =    (RO    -    RU/FLOAT(NGRDFN) 

IF    (     DH    .LT.    0.     )    GO    TO    1 

DELTAS  =  DH*FSTRMN 
DX  =  DH/TAMFMANGL) 
DSTSTR  =  2.*DX 
NPTS  =  1. /DSTSTR  +  1 

IF  (  NPTS  .GT.  400  )  GO  TO  2 


BETA  =  SQRT(FSTRMN*FSTRMN  -  1.) 

RN  =  FLOAT(N) 

RM  =  FLOAT (M) 

RM1  =  FLOAT(Ml) 

FACT1  =  REDFRQ*REDFRQ 

FACT2  =  RN*RN/(FSTRMN*FSTRMN) 
C 

C      COMPUTE  CONSTANTS  TO  BE  USED  IN  THE  COMPUTATIONAL 
C      MOLECULES. 


Al  =  0.25*DELTAS/FSTRMN 
A2  =  0.5*FSTRMN/BETA 
A3  =  REDFRQ*DELTAS*A2 
A4  =  0.5*FSTRMiM 
RETURN 

1  IER  =  1 
WRITE(6,100) 
GO  TO  3 

2  IER  =  I 
WRITE(6,110) 

3  RETURN 

80  FORMAT( «1«  ) 
100  FORMAT*//'  »  ,35X,' ABNORMAL  TERMINATION:  INPUT  RADII  •, 

1    'INCORRECT',//) 
110  FORMATi//8  J  ,35X, 'ABNORMAL  TERMINATION:  MUST  INCREAS', 
1    *E  THE  SIZE  OF  VECTORS  XPT,  CPR,  AND  CPIS,//) 
END 


C  THIS  IS  SUBROUTINE  INITAL 

SUBROUTINE  INITAL 
C 

C      SUBROUTINE  INITAL  INITIALIZES  ALL  FLOW  FIELD 
C      QUANTITIES. 
C 

COMPLEX  DPH1DR, DPHIDX 

COMPLEX  PHI (400) ,DPHISl(400)  ,DPH IS 2(400) 

DIMENSION  X(400) ,R(400) , XPT ( 400 ) ,CPR ( 400 ) , CPU  400 ) 
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COMMON/ BLCK1/  FSTRMN , DELTAS , OX, OH, DSTSTR ,RO  ,R I ,REDFRQ 

COMMON/ BLCK2/  Al , A  2 , A3 , A4 , RN, RM , RM1 , FACT  1 , F ACT  2 

COMMON/ 8 LCK3/  NG  RDFiM  ,  N  PTS  ,  LCOUNT  ,  IHAV  EP  ,  N,  M  ,  Ml ,  I  PR  INT 

C0MMCN/BLCK4/  PH I  * DPHI SI ,OPHIS2 

COMMON/ BLCK5/  X, R  ,  XPT  ,  CPR , CPI 

DATA  PI/3.141593/ 
C 

X( 1)  =  0. 

R ( I )  =  Ru 

XPT(l)  =  0. 
C 

C      ALONG  THE  INITIAL  RIGHT  RUNNING  CHARACTERISTIC  PHI  AND 
C      DPHIS2  ARE  ZERO. 
C 

PHim   =  (o.,o.) 

0PHIS2< 1)    =    (0.,0.) 


AA    =    RM#PI/A4 
DPHISK  1)     =    CMPLX(AA,0.) 
DPHI OR    =    A4-DPHIS1 (1 ) 
DPHIDX    =    A2*DPHIS1(1) 
CPR(l)     =    -2.*REAL(DPHIDX) 
CPK1)     =    -2.*AIMAG(DPHIDX) 

IF    (     IPRINT    .NE.     1     )     RETURN 
WRITE(6,9Q) 
WRITE(6,95J 

WRITE  (6, 100)    XI  li  ,  DPHI  SKI)  ,  DPH  I  S2<  1  )  ,  PH  1  (  1 ) , 
1  R( 1) ,DPHIDR, DPHIDX 

WR1TE(6,110)    CPRCLJtCPICil 
RETURN 
90    FORMAT  (  Il',T10,,X«,T22.,,D(PHI}/D(Sn  •  ,T46, 

1  »D(PHI)/D(S2)  «  ,T70,  «PHIM 

95    FORMATS     •  ,T10,  ■  R«  ,  T23,  •  DC  PHI )  /DR»  ,T47,  ■  0(  PHI  )/DX»  ,  /) 
100    FORMATOH    I    ,  2X,  F8  .4,  4X,  3  (  F  12.  5  ,  F9  .  5  ,3Xi  ,/  •     «,4X,F8.4, 

1  4X,3(F12.5,F9.5,3X) ) 

110    FORMATCO1  ,10X,«  INITIAL    VALUE:        INNER    SURFACE    OF", 

1  ■    OUTER    CYLINDER1, /»     • ,30X, • CP(REAL)     =    »,F8.4,10X, 

2  «CP(  IMAG)    =    ■  ,F8.4,//) 
END 


C  THIS    IS    SUBROUTINE    COMPXR 

C 

SUBROUTINE    COMPXR 
C 

C  SUBROUTINE    COMPXR    COMPUTES    THE     (X,R)    COORDINATES    OF    A 

C  GRID    POINT. 

C 

DIMENSION    X(400) ,R(400) , XPT (400) ,CPR{ 400 ) , CPI ( 400 ) 

COMMON/ BLCK1/    FS TRMN , DELTAS , DX, DH, DSTSTR , RO ,R I , REDFRQ 

COMMON/ BLCK3/    NGRDFN,NPTS, LCOUNT , I  HAVE P , N, M ,M1 , 1  PRINT 

C0MM0N/BLCK5/    X, R , XPT ,CPR, CPI 
C 

I    =    IHAVEP    ♦    1 

J    =    IHAVEP 
C 

C  IF    THIS    IS    THE    FIRST    POINT    ON    A    NEW    RIGHT    RUNNING    MACH 

C  LINE,    GO    TO    1. 


C 


.EQ.    0    )     GO    TO    1 


X(I)     =  X(J)     +    DX 

R(I)    =  R(J)    -DH 

RETURN 

R ( I )    =  RO 

xm  =  xt  i)  *  dst 

RETURN 
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C  THIS    IS    SUBROUTINE    MACHLN 

C 

SUBROUTINE  MACHLN 

C 

C      SUBROUTINE  MACHLN  COMPUTES  THE  VALUES  OF  THE  FLOW 

C      FIELD  QUANTITIES  ALONG  THE  INITIAL  RIGHT-RUNNING  MACH 

C      LINE  WHERE  PHI  AND  DPHIS2  ARE  ZERO. 

C 

COMPLEX  A,C 

COMPLEX  DPHIDR, DPHIDX 

COMPLEX  PHI (400) ,  DPHISH400) ,DPH IS 2(400) 

DIMENSION  X(400)  , R ( 400) , XPT (400 ) ,CPR(400) , CPU  400) 

CCMM0N/BLCK1/  FSTRMN, DELTAS , DX, DH, DSTSTR ,R0 ,R1 , REDFRQ 

COMMON/ BLCK2/  A  1 , A 2, A3 , A4 , RN, RM , RM1, FACT  1, FACT2 

COMMON/ BLCK3/  N3RDFN , NPTS, LCOUNT , IHAVEP , N, M ,M1, I  PRINT 

CQMM0N/BLCK4/  PH  I  , DPH I  SI , DPHI S 2 

C0MM0N/BLCK5/  X, R,XPT , CPR, CPI 

C 

I  =  IHAVEP  +    1 
J  =  IHAVEP 

C 

C      CALCULATE  THE  COEFFICIENTS  ASSOCIATED  WITH  THE  RIGHT- 

C      RUNNING  MACH  LINE. 

C      LINE. 


C 


AA1  =  1.  -  Al/RU) 

A  =  CMPLX( AA1,A3) 

CI  =  I.  +  Al/RU) 

C  =  DPHISUJI*CMPLX(Clt-A3) 

PHI(I)  =  ( 0.,0.) 
DPHIS2( I )  =  (0.,0.) 
DPHISK  I  )  =  C/A 

DPHIDR  =   A4*DPHIS1 (I ) 
DPHIDX  =  A2*DPHIS1(I ) 

IF  (  IPRINT  .NE.  1  )  RETURN 

WRITE (6, 100)  X(  I)  .DPHISKI)  ,DPHIS2(I)  ,PHI<  I), 
1    R( I ) , DPHIDR, DPHIDX 
WRITE(6,111) 
RETURN 
100  F0RMAT(3H  M  , 2X, F 8  .4, 4X,3 ( F12. 5 ,F9. 5 ,3X )  ,/ •  «,4X,F£ 

1    4X,3(F12.5,F9.5,3X) ) 
111  FORMATCO'  ) 
END 


C  THIS  IS  SUBROUTINE  RAD1 

C  THIS  IS  SUBROUTINE  RAD1 
SUBROUTINE  RAD1 

C  SUBROUTINE  RADl  COMPUTES  THE  FLOW  FIELD  QUANTITIES  AT 

C  THE  SURFACE  OF  THE  INNER  RADIUS  WHERE  THE  BOUNDARY 

C  CONDITION  IS   DEPENDENT  ON  N  ,  THE  CIRCUMFERENTIAL 

C  MOOE  NUMBER.   IF   N   IS: 

C  EVEN   —   DPHIDR  =  0,  AND  THUS  DPHIS1  =  DPHIS2. 

C  ODD    —   CP  =  0. 


C 


COMPLEX  A , B , C , T  ANR , DPH I DR , D PH I DX . DEL  T A , s AMM A  , I K 
CCMPL  EX  PH I  (  'i  00  )  t  DPH I  S  1  (  400  )  ,  D  PH  1 1^  I  40  : 
DIMENSION  X(400) , R (400 ) ♦ XPT ( 400 ) ,CPR(400 i ,CPI(400 ) 
COMMON/ BLCK1/  FSTRMN, DELTAS t UX, OH, DSTSTR , RO, RI , REDFRQ 
COMMON/ BLCK 2/  Al , A 2 , A3 , A4, RN,RM , RM1, FACT  1 , FACT2 
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C0MM0N/BLCK3/    NGROFN.NPTS . LCOUNT , IHAVEP ,Nt M ,M1 , I  PR  INT 

C0MM0N/BLCK4/    PHI , DPHI S1,DPHIS2 

COMMON/BLCK5/    X , R , XPT , CPR,CPI 
C 

I    =    IHAVEP    +    1 

J    =    IHAVEP 
C 

C  CALCULATE    THE    COEFFICIENTS    ASSOCIATED    WITH    THE    RIGHT- 

C  RUNNING    MACH    LINE. 

C 

AA1    =    1.    -    Al/R( I ) 

A    =    CMPLX(AAlf A3) 

Bl    =    A1/R(I)    -    0.25*DELTAS*DELTAS*(FACTl    - 
1         FACT2/CRC I)*R(I))) 

B    =    CMPLX.(B1,A3) 

CI    =    1.     *■    Ai/RU) 

C2    =    -A1/R(J)    ♦    0.25*DELTAS*0ELTAS*<FACT1    - 
1         FACT2/(R( I)*R(I))) 

C3   =    (2.*FACT1    -    FACT2/(R( J)*R< J))    - 
1  FACT2/(R(  I)*R(I))  )*DELTAS*0.5 

C  =  DPHIS1 (J)*CMPLX(C1,-A3)  ♦  D PHI S2 ( J) *CMPLX (C2,-A3 ) 
1    +  PHK  J)*CMPLX(C3,0.) 
C 

IF  (  MOD(N,2)  .EQ.  0  )  GO  TO  1 
C 

C      APPLY  THE  BOUNDARY  CONDITION  FOR  ODD  VALUES  OF  N. 
C      CP  =  (IK)PHI  ♦  DPHIDX  =  0.   USE  THIS  B.C.,  THE  FINITE 
C      DIFFERENCE  EQN.  FOR  DPHIDX,  AND  THE  EQN.  FOR  THE  RIGHT 
C      RUNNING  MACHLINE  TO  SOLVE  FOR  DPHIS1  AMD  DPHIS2. 
C 

DELTA  =  A2  *  0.5*IK*DELTAS 

GAMMA    =    IK*(     PHI(J)     *■    0. 5*DELT  AS*DPHI  S2  (  J)     ) 
C 

C  SOLVING    FOR    DPHIS1    AND    DPHIS2 

C 

DPHISKI)    =    (     C    +    B*GAMMA/DELTA     )/(    A    -    B*A2/DELTA    ) 

DPHIS2(I)    =    (     C    -    A*DPHIS1(I)     )/B 
C 

DPHIDX    =    A2*(     DPHISKI)     *    DPHIS2(I)     ) 

DPHIDR    =    A4*(     DPHISKI)    -    0PH1S2(I)     ) 

GO    TO    2 
C 
C  APPLY    THE    BOUNDARY    CONDITION    FOR    EVEN    VALUES    OF       N. 


C 


C 

C  PHI. 


C 


OPHISK  I)    =    C/(A+B) 

DPHIS2(  I)    =    DPHISKI  ) 

DPHIDR    =    (0. ,0.) 

DPHIDX    =    A2*(DPHISKIi     +    DPHIS2(D) 


PHKI)     =    PHKJ)     i-    0.5*(DPHIS2(  J)     *   DPHI  S2(  I  )  )  *DELTAS 

IF    (     IPRINT    .NE.     1     )     RETURN 

WRITE (6, 100)     X(  I) , DPHI  SI ( I)  ,DPHI S2 ( I ) , PH  I  ( I), 
1  R( I) , DPHIDR, DPHIDX 

RETURN 
100    FORMATOH    RI  ,  2X,  F  8  .4,4X,  3(  F  12.  5  ,  F9.5  ,3X  )  t  /  •     «,4X,F3.4, 
1         4X,3(F12.5,F9.5,3X) ) 

END 


C  THIS  IS  SUBROUTINE  RAD2 

C 

SUBROUTINE  RAD2 
C 

C      SUBROUTINE  RAD2  COMPUTES  THE  FLOW  FIELD  QUANTITIES  AT 
C      THE  INNER  SURFACE  OF  THE  OUTER  CYLINDER,  WHERE  TWL 
C      FLOW  TANGENCY  CONDITION  PRESCRIBES  DPHIDR. 
C 
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CCMPLEX  D,EtF,TANR 

CCMPLEX  DPHIDKtDPHIDX 

CCMPLEX  PHI (400)  , DPH I S 1 ( 400  I , DPH  I  S 2 ( 400 ) 

DIMENSION  X(400) ,R(400) ,XPT(400) ,CPR(400) tCPI (400) 

COMMON/ BLCKi/  FSTRMN .DELTA S,DX,DH,DSTSTR,ROtRI ,REDFRQ 

COMMON/ BLCK2/  Al , A 2 , A3 , A4 , RNf KM , RM1 , FAC T 1 T F ACT 2 

C0MM0N/BLCK3/  NGRDFNtNPTSf LCOUNT, IHAVEP, N,M,Ml, I  PRINT 

C0MM0N/BLCK4/  PHI ,DPHI SltD°HIS2 

C0MM0N/BLCK5/  X, R , XPT , CPR, CPI 

DATA  PI/3. 141593/ 
C 

I  =  IHAVEP  4-  1 

J  =  IHAVEP  ♦  2 

K  =  LCOUNT 

ARG  =  RM*PI*X(I) 
C 
C 

C      CALCULATE  THE  COEFFICIENTS  ASSOCIATED  WITH  THE  LEFT 
C      RUNNING  MACH  LINE. 
C 

Dl  =  Al/Rd)  +  0.25*DELTAS*DELTAS*(FACT1  -  FACT2/(R(I) 
1    *R(I))) 

D  =  CMPLX(Dlr~A3) 

El  =  -1.  -  Al/R( I) 

E  =  CMPLX(E1,-A31 

Fi    =    -Al/RU)    -    0.25*DELTAS*DELTAS*{FACTi    -    FACT2/ 
1  (R(I)*R(IJ)I 

F2    =    -1.     *-    A1/R(J) 

F3    =    (2.*FACT1    -    FACT2/(R( I ) *R ( I)J    -    FACT2/(R(J)* 
1  R( J) ) )*DELTAS*0.5 

F    =    DPHISi(J)*CMPLX(Fl,A3)     +    DPH  I  S2(  J)  *CMPLX(  F2,  A3  )    - 
1         PHH  J)*CMPLX(F3,0.) 
C 

C      CALCULATE  THE  FLOW  TANGENCY  CONDITION   -   TANR 
C 
C      BOUNDARY  CONDITIONS  FOR  PANEL  FLUTTER 


C 


C 


TANR1    =    RM*PI*COS(ARG) 
TANR2    =    REDFPvQ*SIN(ARG) 
TANR    =    CMPLX(TANR1 »TANR2i 


DPHIS1U)     =    (     F    +    TANR*E/A4    )/<     D    4-    E    ) 

DPHIS2U)    =    (    F    -    TANR*D/A4    )/(     D    *    E    ) 

DPHIDR    =    A4*(DPHIS1( I)    -    DPHIS2(I)) 

DPHIDX    =    A2*(DPHIS1(I)     +    DPHIS2(I>) 
C 

C  INTEGRATE    ALONG    THE    LEFT    RUNNING    MACH    LINE    TO    FIND    PHI 

C 

PHI(I)     =    PHI(J)     +    0.5*(DPHIS1( J)     *    DPHISKI  >)*DELTAS 
C 
C  CALCULATE    THE    PRESSURE    COEFFICIENT    AT    THIS    GRID    POINT. 


XPT(K)  =  X(I) 

CPR(K)  =  -2.*(-XEDFRQ*AIMAG(PHI( I) )  +  RE AL ( DPH IQX) ] 

CPI(K)  =  -2.*<REDFRQ*REAL(PHI<  I))  ••  AIMAG(  DPHIDX)  ) 


IF  (  IPRINT  .NE.  1  )  RETURN 

WRITE (6, 100)  X(  I ) , DPH  I  SI ( I )  ,DPH I S2( I ) , PH  I( I  )  , 
1    R(I) ,DPHIDR»DPHIDX, TANR 
WRITE(6tllOJ  CPR(K),CPI(K) 
RETURN 
100  F0RMAT(3H  RO, 2X, F8 .4,4X, 3( F 12. 5  ,F9  .  5  ,3X )  ,/  •  •,4X»F8.4, 

1    4X,3(F12.5,F9.5,3X) ) 
110  FORMAT( e0« ,iOX,» INNER  SURFACE  OF  OUTER  CYLINDER',/1  •, 
1    30X, •CP(REAL)  =  » ,F8.4,10X,«CP( IMAG)  =  ',F8.4t//) 
END 
C  THIS  IS  SUBROUTINE  GENFPT 

C 

SUBROUTINE  GENFPTi IER) 
C 

C      SUBROUTINE  GENFPT  COMPUTES  THE  FLOW  FIELD  QUANTITIES 
C      AT  A  GENERAL  FLOW  FIELD  POINT. 
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CCMPLEX  PHIL 

COMPLEX  AfBtCtD, E,F, DENOM 

COMPLEX  DPHIDR, DPHIDX 

COMPLEX  PHI (400)  , DPHI SI (400)  ,DPH  IS  2(400) 

DIMENSION  X(400) tR(400) iXPT(400) tCPR(400) tCPl (400) 

COMMON/ BLCK1/  FS TRMN , DELTAS , DX, DH, DSTST R, RO , RI , REDF RQ 

COMMON/ BLCK2/  A  1 , A2 , A3 , A4 ,R  N , RM , RM1 , F ACT 1 , F ACT2 

COMMON/ BLCK3/  NGRDFN , NPTS , LC OUNT , IHAVEP , N, M , Ml , I  PR  INT 

COMMON/ BLCK4/  PHI,DPHI S1,DPHIS2 

C0MM0N/BLCK5/  X , R , XPT , CPR , CP I 
C 

I  =  IHAVEP  +  1 

J  =  IHAVEP 

K  =  IHAVEP  +  2 
C 
C      COMPUTE  THE  COEFFICIENTS  OF  DPHIS1  AND  DPHIS2. 

FACT3  =  FACT2/(R(I )*R< I)) 

FACT4    =    FACT2/(R(K)*R(K) ) 

FACT5    =    FA'CT2/(R(  J)*R(  J)) 
C 

C  A,B    AND    C    ARE    ASSOCIATED    WITH    THE    RIGHT    RUNNING    MACH 

C  LINE,     WHILE    D,E    AND    F    ARE    ASSOCIATED    WITH    THE    LEFT. 

AA1    =    1.    -    Al/R( I) 

A    =    CMPLX(AA1, A3) 

Bl    =    -AA1    +    1.     -    0.25*DELTAS*DELTAS*(FACT1    -    FACT3) 

B    =    CMPLX(B1,A3) 

CI   =    1.    «■    A1/R(J) 

C2    =    -CI    *■    1.    +    0.25*DELTAS*DELTAS*(FACT1    -    FACT3) 

C3   =    (2.*FACT1    -    FACT3    -    FACT5 ) *DELTAS*0 .5 

C    =    DPHIS1(J)*CMPLX(C1,-A3)     *    DPHI S2< J) *CMPLX(C2 ,-A3) 
1  +    PHIU)*CMPLX(C3,0.) 

C 

Dl    =    A1/R(I)    <-    0.25*DELTAS*DELTAS*(FACT1    -    FACT3) 

D    =    CMPLX(Di,-A3) 

El    =    -1.    -    Al/RCI) 

E    =    CMPLX( El, -A3) 

Fl    =    -A1/R(K)    -    0.25*DELTAS*DELTAS*(FACT1    -    FACT3) 

F2    =    -1.    +    A1/R(K) 

F3    =    (2.*FACT1    -    FACT3    -    F ACT4) *DELTAS*0 . 5 

F    =    DPHISHK)*CMPLX(F1,A3)     +    D PHIS2IK) *CMPLX(F2, A3)    - 
1         PHI(K)*CMPLX(F3,0.) 

C  USE    CRAMERS    RULE    TO    SOLVE    FOR    DPHIS1    AND    DPHIS2 

C  PROVIDED    THAT    THE    DETERMINANT    OF    COEFFICIENTS    IS    NON- 

C  SINGULAR. 

C 

DENOM    =    A*E   -    D*B 

IF    (     CABS(     DENOM    )     .LT.    l.E-05     )     GO    TO    1 

DPHISKI)    =    (C*E    -    B*F)/DENOM 

DPHIS2CI)     =    (A*F    -    C*D    J/DENOM 

DPHIDR    =    A4*(DPHIS1(I)    -    DPHIS2(D) 

DPHIDX    =    A2*(DPHIS1(I)     +    DPHIS2(I)) 
C 

C  INTEGRATE    ALONG    THE    RIGHT    RUNNING    MACH    LINE    TO    FIND 

C  PHI. 


C 


PHI(I)  =  PHI(J)  *  0.5*(DPHIS2(  I  )  +-  DPHI  S2  (  J  i  )  ^DELTAS 
PHIL  =  PH-KKJ  +    0.5*(DPHIS1(K)  ♦  DPHIS 1  (  I)  }  ^DELTAS 

IF  (  IPRINT  .NE.  1  )  RETURN 

WRITE (6, 100)  X( I ) , DPHI  SI ( I ) , DPHI S2( I ) , PHI ( I), 

1    R( I) , DPHIDR, DPHIDX, PHIL 
WRITE(6,111) 
RETURN 

1  IER  =  1 

WRITE( 6,110)  DENOM 
RETURN 
100  FGRMATC3H  G  , 2X, F 8 .4 , 4X, 3 ( F 12. 5 , F9 .5 ,3X ) , / •  ',4X,F8.4, 
1    4X,3(F12.5,F9.5,3X) ) 
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MAT(  •()»  ,///,•  ■ »35X, "ABNORMAL  TERMINATION:  •, 
•DETERMINANT  OF  COEFFICIENTS  IS  SINGULAR  •  ,/■()•,  ! 
•DENOMINATOR  =  «,E12.5) 


110  FORMAT( «0»  ,///,' 

1  'DETERMINANT  OF  COEFFICIENTS  IS  SINGULAR  •  ,/■()•,  50X, 

2  •OENOMIN "~ 

111  FORMAT(«0»  ) 
END 


C  THIS  IS  SUBROUTINE  QSF 

C         SUBROUTINE  QSF 

SUBROUTINE  QSF ( H , Y , Z , NDIM ) 
C  CALL  QSF  (H,Y,Z,NDIM) 

C         DESCRIPTION  OF  PARAMETERS 

C  H         THE  INCREMENT  OF  ARGUMENT  VALUES. 

C  Y         THE  INPUT  VECTOR  OF  FUNCTION  VALUES.. 

C  Z         THE  RESULTING  VECTOR  OF  INTEGRAL 

C  VALUES.   Z   MAY  BE  IDENTICAL  WITH   Y. 

C  IDENTICAL  WITH  Y. 

C  NDIM    -  THE  DIMENSION  OF  VECTORS  Y  AND  Z. 

C 

C         REMARKS 

C  NO  ACTION  IN  CASE  NDIM  LESS  THAN  3. 

C 
C 
C 

DIMENSION  Y(i) ,Z(1) 
C 

HT=.3333333*H 

IF(NDIM-5)7,8,i 
C 

C      NDIM  IS  GREATER  THAN  5.  PREPARATIONS  OF  INTEGRATION 
C      LOOP 
C 

1  SUM1=Y(2J*Y(2) 
SUM1=SUM1+SUM1 
SUM1=HT*(Y(1  )+SUMH-Y(3)) 
AUX1=Y(4)+Y(4) 
AUX1=AUXH-AUX1 

AUX1  =  SUMH-HT*(Y(3)+AUXH-Y(5)  ) 

AUX2=HT*(Y(1)  +  3.875*(Y(2)+Y(5)H-2.625*{Y(3)*Y(4))  + 
1    Y(6)) 
SUM2=Y( 5)*Y(5) 
SUM2=SUM2+SUM2 

SUM2=AUX2-HT* ( Y ( 4) +SUM2+Y{ 6 ) ) 
Z(1)=0. 
AUX=Y(3)*Y(3) 
AUX=AUX+AUX 

Z(2)=SUM2-HT*(Y(  2)+AUX+Y(4)  ) 
Z(3)=SUM1 
Z(4)=SUM2 
IF(NDIM-6)5,5,2 

C      INTEGRATION  LOOP 

2  DO  4  1  =  7, NDIM,  2 
SUM1=AUX1 
SUM2=AUX2 

AUX1  =  Y(  I-l)fy(I-l) 

AUX1=AUX1.+  AUX1 

AUX1  =  SUMH-HT*(Y<  1-2) + AUXi+Y < I )  ) 

Z(I-2)=SUM1 

IF(I-NDIM)3,6, 6 

3  AUX2=Y(  I  )*Y(  I) 
AUX2  =  AUX2.+  AUX2 
AUX2=SUM2*HT*(Y(  I- 1 ) +AUX2 *Y <  IU)J 

4  Z(I-1)=SUM2 

5  Z(NDIM-1)=AUX1 
Z(NDIM)=AUX2 
RETURN 

6  Z(NU1M-1)=SUM2 
Z(NDIM)=AUX1 
RETURN 
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C  END    OF     INTEGRATION    LOOP 

C 

7  IFCNDIM-3) 12,11,8 
C 

C      NDIM  IS  EQUAL  TO  4  OR  5 

8  SUM  2=1.  12  5*HT*(  YU  )  +  Y(  2) +Y<  2)+Y<  2) +Y(3)+Y<  3  )  +  Y(  3)  * 
1    Y(4) ) 

SUM1=Y( 2)+Y(2) 

SUM1  =  SUM1«-SUM1 

SUMl=HT*(Ym+SUMH-Y(3)  ) 

Z(U=0. 

AUXl=Y<  3)+Y<3) 

AUXi  =  AUXl.+  AUXl 

Z(2)  =  SUM2-HT*(Y(2H-AUXH-Y<4)) 

IF(NDIM-5)10,9,9 

9  AUX1=Y(4)+Y<4) 
AUXl  =  AUXH-AUXi 
Z(5)=SUM1+HT*(Y(  3)*AUXH-Y(5)) 

10  Z(3)=SUM1 
Z(4)=SUM2 
RETURN 

C 

C  NDIM    IS    EQUAL    TO    3 

11  SUM1=HT*(1.2  5*Y< 1 ) *Y( 2 ) +Y( 2)-. 25*Y<3) ) 
SUM2=Y(2)+Y(2) 

SUM2=SUM2+SUM2 

Z(3)=HT*<Y<1)+SUM2+Y(3)) 

Z(1J=0. 

Z(2)=SUM1 

12  RETURN 
END 
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